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ABSTRACT 


An  experimental  test  of  a  structure  results  in  a  database  of  frequency  response 
functions  from  which  natural  frequencies  and  mode  shapes  are  identified  for  the  given 
boundary  conditions  of  the  test.  The  natural  frequencies  of  a  system  under  a  variety  of 
boundary  conditions  can  be  identified  by  applying  artificial  boundary  conditions  (ABCs) 
at  measurement  locations  and  obtaining  the  frequency  response  function  for  the 
unrestrained  degrees  of  freedom.  These  frequencies  are  found  without  any  physical 
alterations  of  the  test  boundary  conditions.  Use  of  the  ABCs  in  sensitivity-based  model 
updating  and  damage  detection  in  conjunction  with  baseline  data  provides  improved  error 
localization  and  results  in  a  more  accurate  finite  element  model. 


vi 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . 1 

II.  THEORY . 5 

A.  OMITTED-COORDINATE  SET . 5 

B.  EXACT  DYNAMIC  REDUCTION  AND  FRF  MATRICES . 6 

C.  DRIVING  POINT  FREQUENCY  RESPONSE  FUNCTION . 7 

III.  ABC  CONFIGURATION  FREQUENCIES . 1 1 

A.  DRIVING  POINT  ANTIRESONANCES  DEFINE  ABC 

FREQUENCIES  FOR  A  SINGLE  A-SET  COORDINATE . 11 

B.  ABC  CALCULATIONS  FOR  A  2  DOF  SYSTEM . 12 

C.  ABC  FREQUENCY  CALCULATIONS  FOR  FREE-FREE  BEAM  14 

IV.  ABC  CONFIGURATION  IN  SENSITIVITY-BASED  UPDATING . 19 

A.  SENSITIVITY  MATRIX  DEFINED . 19 

B.  IMPROVED  CONDITIONING  USING  ABC  SENSITIVITIES . 21 

C.  COMPUTER  SIMULATED  DAMAGE  DETECTION  USING  ABC 

SYSTEM  SENSITIVITIES . 21 

1.  Damage  Detection:  Example  1 . 22 

2.  Damage  Detection:  Example  2 . 23 

3.  Damage  Detection:  Example  3 . 25 

4.  Damage  Detection:  Example  4 . 25 

5.  Damage  Detection:  Example  5 . 28 

6.  Damage  Detection:  Example  6 . 29 

V.  EXPERIMENTAL  APPLICATION . 31 

vii 


A.  EQUIPMENT  SETUP . 31 

B.  DATA  COLLECTION . 33 

C.  O-SET  FREQUENCIES  USING  EXPERIMENTAL  DATA . 33 

D.  DAMAGE  DETECTION  USING  EXPERIMENTAL  DATA . 43 

VI.  CONCLUSIONS  AND  RECOMMENDATIONS . 47 

A.  CONCLUSIONS . 47 

B.  RECOMMENDATIONS . 48 

APPENDIX  MATLAB  CODE . 49 

A.  MAIN  PROGRAMS . 49 

B.  SUPPORT  PROGRAMS . 64 

C.  FUNCTIONS  CALLED  BY  MAIN  &  SUPPORT  PROGRAMS . 73 

LIST  OF  REFERENCES . 77 

INITIAL  DISTRIBUTION  LIST . 79 


viii 


LIST  OF  FIGURES 


Figure  3.1  2  DOF  System . 1 1 

Figure  3.2  2  DOF  System  with  DOF  #1  restrained  to  ground . 12 

Figure  3.3  Hn(Q)  versus  Q . — •  13 

Figure  3.4  H],  (Q)  versus  Q . 14 

Figure  3.5  Free-free  beam . 15 

Figure  3.6  Driving  Point  FRF  Hn(Q)  for  Free-Free  beam . 15 

Figure  3.7  H]]1  (Q)  for  free-free  beam  (a-set  DOFs:  1,5,9,13,17,21) . 16 

Figure  3.8  ABC  Configuration . 17 

Figure  3.9  H^1  (Q)  for  free-free  beam  (a-set  DOFs:  1, 17) . 18 

Figure  3.10  ABC  Configuration  -  a-set:  [1  17] . 18 

Figure  4.1  Damage  Detection  -  Example  1 . 23 

Figure  4.2  Damage  Detection  -  Example  2 . 24 

Figure  4.3  Damage  Detection  -  Example  3 . 26 

Figure  4.4  Damage  Detection  -  Example  4 . 27 

Figure  4.5  Damage  Detection  -  Example  5 . 28 

Figure  4.6  Damage  Detection  -  Example  6 . 29 

Figure  5.1  Laboratory  Setup . 32 

Figure  5.2  Hi  i(Q)  for  free-free  beam . 34 

Figure  5.3  ABC  Configuration:  Node  1  restrained  to  ground . 35 

Figure  5.4  H;]  (£2)  for  Node  1  restrained  to  ground . 36 

Figure  5.5  ABC  Configuration:  Nodes  1  &  2  restrained  to  ground . 37 

Figure  5.6  H~]  (Q)  for  Nodes  1,2  restrained  to  ground . 37 

Figure  5.7  ABC  Configuration:  Nodes  1,4  restrained  to  ground . 38 

Figure  5.8  H (Q)  for  Nodes  1,4  restrained  to  ground . 38 

Figure  5.9  ABC  Configuration:  Nodes  1,2,3  restrained  to  ground . 39 

Figure  5.10  H~](Q)  for  Nodes  1,2,3  restrained  to  ground . 39 

Figure  5.1 1  ABC  Configuration:  Nodes  1,2, 3, 4  restrained  to  ground . 40 


IX 


Figure  5.12  H;](Q)  for  Nodes  1,2, 3, 4  restrained  to  ground . 40 

Figure  5.13  ABC  Configuration:  Nodes  1,2,3, 4, 5  restrained  to  ground . 41 

Figure  5.14  H'](Q)  for  Nodes  1,2, 3, 4, 5  restrained  to  ground . 41 

Figure  5.15  ABC  Configuration:  Nodes  1,3,5  restrained  to  ground . 42 

Figure  5.16  H;]  (Q)  for  Nodes  1,3,5  restrained  to  ground . 42 

Figure  5.17  Error  localization  using  FE  frequencies . 43 

Figure  5.18  Error  localization  using  test  frequencies . 44 


LIST  OF  TABLES 


Table  4. 1  Example  1 :  FE  and  Test  Frequencies  and  Error . 22 

Table  4.2  Example  2:  FE  and  Test  Frequencies  and  Error . 24 

Table  4.3  Example  3:  FE  and  Test  Frequencies  and  Error . 25 

Table  4.4  Example  4:  FE  and  Test  Frequencies  and  Error . 27 

Table  5.1  FE  &  Test  Frequencies  for  free-free  beam . 34 

Table  5.2  FE  &  Test  Frequencies  for  Node  1  restrained  to  ground . 35 

Table  5.3  FE  &  Test  Frequencies  for  Nodes  1,2  restrained  to  ground . 37 

Table  5.4  FE  &  Test  Frequencies  for  Nodes  1,4  restrained  to  ground . 38 

Table  5.5  FE  &  Test  Frequencies  for  Nodes  1,2,3  restrained  to  ground . 39 

Table  5.6  FE  &  Test  Frequencies  for  Nodes  1,2, 3, 4  restrained  to  ground . 40 

Table  5.7  FE  &  Test  Frequencies  for  Nodes  1,2, 3, 4,5  restrained  to  ground . 41 

Table  5.8  FE  &  Test  Frequencies  for  Nodes  1,3,5  restrained  to  ground . 42 


XI 


xii 


ACKNOWLEDGMENTS 


I  would  like  to  express  my  sincerest  appreciation  to  Professor  Joshua  H.  Gordis 
for  his  dedicated  support  and  undying  patience  throughout  this  endeavor. 

I  would  also  like  to  thank  Frank  Franzen  for  his  craftsmanship  and  Tom  Christian 
for  never  giving  up  on  the  installation  of  I-DEAS  and  making  the  laboratory  experiment  a 
success. 


xm 


I.  INTRODUCTION 


In  general,  a  finite  element  (FE)  model  is  defined  by  a  large  number  of  physical 
parameters,  but  a  modal  test  of  a  structure  results  in  a  small  number  of  modal  parameters 
which  are  used  to  modify  the  model  parameters.  It  is  always  crucial,  then,  to  improve  or 
"update"  a  finite  element  model  in  order  to  produce  a  reliable,  confident  prediction  of  the 
structural  response  of  a  system.  Adjusting  the  parameters  that  define  the  model  can 
reduce  inaccuracies  in  a  FE  model.  These  parameters  include,  but  are  not  limited  to 
dimensional  properties,  moduli  of  elasticity,  and  density.  If  the  parameters  in  the  model 
can  be  adjusted  for  closer  agreement  with  the  measured  parameters  (identified  in  a  modal 
test),  then  the  model  will  provide  a  more  accurate  representation  of  dynamic  response.  It 
is  important  to  remember  that  no  FE  model  can  exactly  predict  the  dynamic  response  of  a 
real  structure,  but  can  still  provide  reasonable  and  reliable  results.  This  is  the  goal  of 
model  updating. 

One  problem  that  has  been  historically  associated  with  finite  element  model 
verification  has  been  in  measuring  enough  information.  Procedures  have  been  developed 
for  measuring  larger  experimental  databases.  One  such  method  that  has  been  developed 
is  known  as  "Perturbed  Boundary  Condition"  testing.  [Ref.  4]  In  this  procedure,  selected 
boundary  conditions  of  a  system  under  test  are  drastically  changed  and  the  system  is 
retested.  This  method  is  effective,  however  the  procedure  requires  that  physical 
modifications  be  made  to  the  structure  and  additional  tests  are  necessary  for  each 
modification  made.  The  result  is  increased  time,  increased  work  and  most  importantly  in 
today's  world,  increased  cost. 
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A  more  efficient  method  is  to  identify  additional  and  distinct  mode  frequencies 
from  the  same  modal  test  that  was  performed  to  identify  the  mode  frequencies  of  the 
standard  system,  without  the  need  for  physical  modification  of  the  structure.  The 
additional  frequencies  and  mode  shapes  correspond  exactly  to  the  mode  frequencies 
found  when  combinations  of  measured  coordinates  are  constrained  to  ground.  These 
additional  frequencies  are  therefore  associated  with  different  boundary  conditions  for  the 
structure  and  no  physical  change  in  the  boundary  conditions  has  been  made.  [Ref.  1] 
Because  the  boundary  conditions  are  not  actually  applied,  they  are  termed 
"artificial  boundary  conditons",  or  "ABCs."  For  a  single  degree  of  freedom  system  the 
driving  point  antiresonance  frequencies  of  a  frequency  response  function  (FRF) 
correspond  to  the  frequencies  of  the  structure  with  the  driving  point  degree-of-ffeedom 
restrained  to  ground.  The  ABCs  are  imposed  on  the  FE  model  and  correspond  to  ideal 
constraints.  With  only  one  experimental  database,  this  scheme  yields  a  separate  FE 
model  for  each  ABC  configuration.  With  the  exception  of  the  imposed  boundary 
conditions,  each  model  is  identical  and  can  be  used  to  generate  sensitivity  data.  In  a  set 
of  spatially  incomplete  frequency  response  function  data,  the  ABCs  are  the  boundary 
conditions  that  define  an  "omitted  coordinate  system  (OCS),"  or  "o-set."  In  performing  a 
vibration  test,  a  choice  is  made  as  to  the  set  of  coordinates  to  instrument  with  response 
transducers.  This  is  the  "analysis  set,"  or  the  "a-set."  With  the  choice  of  an  a-set,  the  o- 
set  is  defined  as  the  complimentary  set  of  coordinates.  For  an  FE  model,  the  o-set  is  of 
finite  dimension,  but  with  respect  to  the  test  system,  this  set  is  of  infinite  dimension.  The 
natural  mode  frequencies  of  the  o-set  systems  appear  as  peaks  in  the  impedance  spectra. 
[Ref.  2] 
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ABC  configuration  frequencies  are  available  from  any  set  of  test  data  due  to  the 
fact  that  a  spatially  incomplete  FRF  matrix  is  identically  equal  to  the  FRF  matrix  that  is 
calculated  from  the  exact  dynamic  reduction.  [Ref.  1] 

In  addition  to  providing  a  larger  number  of  frequencies  for  a  system  using  one 
experimental  database,  the  ABC  also  provides  a  means  to  eliminate  or  greatly  reduce  ill- 
conditioning  in  sensitivity  equations.  With  the  system  artificially  restrained  at  various 
measured  coordinates,  the  columns  of  the  sensitivity  matrix  which  are  found  from  the 
ABC  configuration  are  linearly  independent  of  the  columns  calculated  from  the  baseline 
configuration.  This  reduces  the  difficulties  in  determining  the  parameters  that  are  in 
error.  Most  importantly,  in  order  to  more  accurately  update  the  FE  model,  this 
configuration  improves  error  localization. 
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II.  THEORY 


A.  OMITTED-COORDINATE  SET 

The  equation  for  the  steady-state  harmonic  response  for  a  full  order  model 

(without  damping)  is 


]Xa 

kao' 

-Q2 

maa 

3 

o 

|Xal  _|fa] 

_ka 
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.moa 

mooJ 

,lXoJ  lfoJ 

where  Q  is  the  frequency  of  harmonic  excitation  and  f  is  a  vector  of  excitations  consistent 
with  x  and  k  and  m  are  the  stiffness  and  mass  matrices,  respectively.  The  subscript  "a" 
refers  to  the  "a-set"  and  the  subscript  "o"  refers  to  the  "o-set".  Equation  2.1  is  also  known 
as  the  impedance  model  and  can  be  redefined  as: 


(2.2) 


where  Z(Q)  is  the  system  impedance  matrix  evaluated  at  the  forcing  frequency,  Q. 
Assuming  that  {f0}  =  {0}  (static  constraint  -  no  excitation  acting  on  the  omitted 
coordinates),  the  exact  relationship  between  the  o-set  and  a-set  coordinates  is  as  follows: 

{x0}  =  [l-O2k;>00]"'[-  k^koa  +  Q2k;>oa]{xa}  (2.3a) 

or 

KH-zJKIM  '  <2-3b> 


The  origin  of  the  o-set  system  results  from  Eq.  (2.3a).  By  definition,  the  bracketed 
inverse  term  is 
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where  Det[  •  ]  indicates  the  determinant  and  Adj[  •  ]  indicates  the  adjoint  matrix.  From 
Eq.  (2.4),  it  is  clear  that  the  bracketed  inverse  term  does  not  exist  at  those  frequencies  £2 
which  satisfy 

Det[l-Q2k^moo]  =  0  (2.5) 


Therefore  the  relationship  between  the  a-set  and  o-set  does  not  exist  at  those  frequencies. 
The  frequencies  that  satisfy  Eq.  (2.5)  are  the  eigenvalues  of  the  system  defined  by  [k00] 
and  [iHoo],  the  o-set  system.  This  system  is  obtained  by  fully  constraining  all  coordinates 
in  the  a-set  to  ground.  The  use  of  spatially  incomplete  FRF  data  in  identification  of  the  o- 
set  mode  frequencies  implies  the  use  of  dynamic  reduction  applied  to  the  structure  under 
test.  [Ref.  2] 

B.  EXACT  DYNAMIC  REDUCTION  AND  FRF  MATRICES 


By  equipping  a  structure  with  a  finite  number  of  response  transducers,  a  reduced 
order  model  is  defined,  where  the  impedance  of  the  reduced  order  model  is  not  linearly 
dependent  on  the  impedance  of  the  full  order  model.  The  "full-order",  exact  FRF  model 
of  a  structure  is  a  FRF  matrix  of  infinite  dimension, 


Hx 


Hao 

Hoa 

H00 

(2.6) 


where  the  number  of  coordinates  in  the  o-set  is  infinite.  The  FRF  matrix  measured  in  a 
test  is  a  matrix  partition  extracted  from  the  infinite  dimension  matrix, 

Hx  =  Haa  (2.7) 
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where  the  overbar  notation  indicates  a  reduced  model.  The  matrix  partition  defined  by 
Eq.  (2.7)  represents  a  structural  dynamic  model  that  has  been  reduced  using  exact 
dynamic  reduction.  From  the  partitioned  inverse  relation  of  the  FRF  matrix  to  its 
associated  impedance,  ZH  =  I.  [Ref.  2] 

H„=(z„-z„z;xr  (2-8) 

The  presence  of  the  o-set  system  dynamics  in  the  term  Z"’  is  seen  in  Eq.  (2.8).  The 
formula  for  the  matrix  inverse  seen  in  Eq.  (2.4)  is  applicable.  Since  every  element  in  Z^ 
is  singular  at  the  natural  frequencies  of  the  o-set,  Eq.  (2.8)  shows  that  elements  of  H'j 
will  be  singular  at  the  o-set  natural  frequencies  [Ref.  1].  Since  the  measured  FRF  matrix 
implicitly  defines  a  dynamically  reduced  impedance  model,  an  order  n  FRF  matrix  is 
calculated  from  the  FE  model,  and  the  partition  of  this  matrix  which  corresponds  to  the 
coordinates  measured  in  a  test  is  extracted  (the  a-set).  The  reduced  model  retains  the 
entire  modal  content  of  the  original  model. 

C.  DRIVING  POINT  FREQUENCY  RESPONSE  FUNCTION 

Starting  again  with  the  equation  of  steady  state  forced  response  for  a  linear 
structural  dynamic  system,  Eq.  (2.2),  rewritten  (with  damping)  as 

[k  -  Q2m  -  jQC]{x}  =  (f }  (2.9) 

where  [k  -  Q2m  -  j£2C]  =  [Z]  and  C  is  the  damping  matrix.  Transferring  from  physical 
to  modal  coordinates  let 

{x}=  [<D]{q}  (2.10) 

Rewriting  Eq.  (2.9), 
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[Z][0]{q}  =  {f} 


(2.11) 


Premultiply  by  [0]T , 

[0]T[Z][<b]{q}  =  [0]T{f}  (2.12a) 

[0Tk0-Q20Tm0  +  jQOTCO]{q}  =  {3}  (2.12b) 


Using  orthogonality  and  assuming  proportional  damping, 

[co2  - Q2  +  2j^Qco]{q}  =  {3}  (2.13) 

where  to,  is  the  natural  frequency  of  the  ith  mode,  £  is  the  damping  ratio  and 
[co2  -Q2  +  2j£Qco]  is  a  diagonal  matrix.  This  matrix  is  inverted  to  find  the  modal 
frequency  response: 


{q}  = 


CD2  -Q  +  2}C)Q(£>i] 


{3} 


(2.14) 


Transforming  back  into  physical  coordinates  by  premultiplying  by  [O]  and  using 

{3}  =  [0]T{f},  the  modal  decomposition  of  the  FRF  matrix  in  physical  coordinates  is  as 

follows: 


W=[o] 


|_(o2  -  Or  +2j^Qc0j  J 


[o]T{f} 


(2.15) 


where  [H(Q)]  =  [O] 


1 


to2  -Q2+2j;Qcoj 


[Of .  [H(Q)]  can  be  written: 


NDOF 


[H(Q)]=  Z 


MM1 


t  co2  -Q  +2j^Qcoj 


(2.16a) 
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or  any  element, 


NDOF 

Hij(O)  -  £ 

P-1  ®i 


{♦!•}{♦?} 

-Q2+2j^Qcoi 


(2.16b) 


This  driving  point  FRF  will  be  used  in  subsequent  chapters  in  defining  the  ABC  natural 
frequencies. 
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III.  ABC  CONFIGURATION  FREQUENCIES 


A.  DRIVING  POINT  ANTIRESONANCES  DEFINE  ABC  FREQUENCIES 
FOR  A  SINGLE  A-SET  COORDINATE 

The  following  example  using  a  basic  2  degree-of-freedom  (DOF)  system,  shows 
that  driving  point  antiresonances  correspond  to  the  natural  frequencies  of  the  structure 
with  the  driving  point  DOF  constrained  to  ground. 


mi  ki 


m2  k2 


Figure  3.1  2  DOF  system 


From  Eq.  (2.16b),  the  undamped  driving  point  FRF  is  given  by 


NDOF 


(3.1) 


where  (j),  is  a  mass  normalized  mode  shape  element,  ©r  is  the  rth  natural  frequency,  and  Q 


is  the  forcing  frequency.  The  frequency  of  the  anti-resonance  of  Hn(Q)  is  given  by 


2  _  Ru©2  +  R-n^i 

^ anti -res  r>  1  t>  2 

Kn  +  Kn 


(3.2) 


where  the  modal  residue  is  given  by  R-j  =  (J)^  - .  It  can  be  demonstrated  that  this 

frequency  is  identical  to  the  natural  frequency  of  the  system  in  Figure  3.1,  with  the 
driving  point  DOF  constrained  to  ground  as  in  Figure  3.2. 
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Figure  3.2  2  DOF  system  with  DOF  #1  restrained  to  ground. 

Using  a  simplified  numerical  example,  if  the  values  for  Figure  3.1  are  set  as  mi=m2=l  .0 
and  ki=k2=1.0,  the  frequency  of  the  single  antiresonance  is  Qanti-res =  42  rad/sec,  which 
is  identically  equal  to  the  single  natural  frequency  of  the  ABC  system  seen  in  Figure  3.2, 

which  is  co  =  42  rad/sec. 

B.  ABC  CALCULATIONS  FOR  A  2  DOF  SYSTEM 

The  ABC  frequencies  are  calculated  using  Eq.  (2.8).  A  spatially  incomplete  FRF 
matrix  is  generated  and  inverted  at  each  frequency.  In  a  plot  of  the  elements  of  the 
resulting  impedance  matrix  versus  frequency,  the  singular  frequencies  are  also  the  ABC 
system  frequencies.  To  show  an  example,  the  2  DOF  system  shown  in  Figure  3.1  is  used 
to  calculate  the  driving  point  FRF  Hi  i(O).  The  system  contains  two  modes  and  the 
frequencies  of  these  modes  are  coi  =  0.618  rad/sec  (fj  =  0.0984  Hz)  and  co2  =  1  -618 
rad/sec  (f2  =  0.2575  Hz),  respectively.  A  plot  of  HU(Q)  is  shown  in  Figure  3.3  and  these 
modes  can  be  seen  on  the  plot.  In  addition  to  these  two  modes,  it  can  also  be  seen  that 
the  anti-resonance  frequency,  as  computed  in  Eq.  (3.2),  is  Qami-res =  42  =  1.414  rad/sec. 
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Figure  3.3  Hn(Q)  versus  Q 

Using  Eq.  (2.8),  the  ABC  system  frequency  can  be  found  by  calculating  H']  (Q)  where 
the  a-set  is  defined  as  DOF  #1,  and  (O)  is  the  scalar  inverse  of  Hn(Q).  H;](0)  is 
plotted  versus  Q  in  Figure  3.4.  The  single  singular  frequency  corresponds  to  the  anti¬ 
resonance  frequency  for  this  case  (QWi-res  =1.414  rad/sec).  This  singular  frequency  is 
also  the  natural  frequency  of  the  ABC  system  that  was  obtained  by  constraining  DOF  #1 
to  ground  (as  in  Figure  3.2).  Most  importantly,  this  simple  example  demonstrates  the 
power  and  usefulness  of  Eq.  (2.8)  and  shows  that  driving  point  antiresonances  correspond 
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Figure  3.4  H,,1  (Q)  versus  Q 
(Natural  Frequency  of  ABC  system) 

directly  to  the  natural  frequencies  of  the  system  with  the  driving  point  coordinate 
constrained  to  ground.  [Ref.  1] 

C.  ABC  FREQUENCY  CALCULATIONS  FOR  A  FREE-FREE  BEAM 

To  further  demonstrate  the  usefulness  of  the  ABC  technique,  a  free-free  beam  is 
examined.  The  model  consists  of  10  beam  elements  and  each  element  contains  two 
nodes.  The  nodes  have  2  degrees-of-freedom,  translational  and  rotational,  giving  the 
system  a  total  of  22  DOF.  The  odd-numbered  DOF  are  the  translational  degrees  and  the 
even-numbered  DOF  are  the  rotational  degrees.  The  properties  of  the  beam  are  as 
follows:  length  =  60  inches,  El  =  5.5x10s  lbf-in2,  p  =  0.283  lbf/in2,  and  cross-sectional 
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area  =  0.75  in2.  Response  transducers  (accelerometers)  are  placed  at  every  other  node,  at 
the  translational  DOFs  #1,5,9,13,17,  and  21.  This  model  is  shown  in  Figure  3.5 

----------- 

J  _ 

Figure  3.5  Free-free  beam  -  Transducer:  [_] 

The  set  of  coordinates  chosen  to  be  instrumented  with  response  transducers  is  the  a-set,  in 
this  example  defined  as  [1  5913  1721].  It  can  be  assumed  by  the  given  setup,  that 
excitation  is  applied  at  each  of  the  a-set  DOFs,  resulting  in  a  6  by  6  matrix.  The 
impedance  matrix,  H;]  (Q),  the  scalar  inverse  of  Hn(Q),  is  calculated  at  excitation 
frequencies  from  0-800  Hz.  Figure  3.6  shows  the  driving  point  FRF  for  the  system. 


Frequency  (Hz) 


Figure  3.6  Driving  Point  FRF  Hn(Q)  for  Free-Free  Beam 
(a-set  DOFs:  1,5,9,13,17,21) 
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The  ABC  system  frequencies  are  found  by  calculating  the  impedance  matrix,  H;](Q). 
By  plotting  the  1,1  element  (H  7/  (P))  of  the  matrix,  the  ABC  system  frequencies  can  be 
seen  (Figure  3.7). 


Figure  3.7  H (Q)  for  Free-Free  Beam 
(a-setDOFs:  1,5,9,13,17,21) 

These  ABC  frequencies  correspond  exactly  to  the  natural  frequencies  of  the  system  with 
all  of  the  measured  coordinates  (a-set)  constrained  to  ground.  This  configuration  is 
shown  in  Figure  3.8.  Confirming  the  results,  the  natural  frequencies  of  a  beam  with 
pin  restraints  in  the  locations  shown  in  Figure  3.8  are  calculated  to  be:  346.70,  385.13, 
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Figure  3.8  ABC  Configuration  -  a-set:  [1  5913  1721] 
(Measured  coordinates  restrained  to  ground) 


482.73,  610.23,  and  735.1 1  Hz.  These  frequencies  correspond  exactly  to  the  peaks 
shown  in  Figure  3.7. 

A  common  vibrational  test  involves  two  shakers  placed  to  excite  all  modes  in  the 
bandwidth  of  interest.  To  simulate  this,  the  next  ABC  configuration  set  will  be  DOFs  1 
and  17,  defined  as  the  translations  at  nodes  1  and  9.  Performing  the  same  calculations  as 
above,  the  impedance  matrix,  H  ~la  (Q),  is  obtained.  Looking  at  the  1 , 1  element  of  this 

matrix,  an  additional  set  of  ABC  frequencies  is  found.  A  plot  of  H  (Q)  over  the  range 
0-800  Hz  is  shown  in  Figure  3.9.  The  equivalent  system  is  a  beam  with  pin  restraints  at 
node  1  and  node  9.  The  natural  frequencies  of  this  system  are  calculated  to  be  20.45, 
63.57, 112.30, 215.05, 364.61,  544.16,  and  669.70  Hz.  These  frequencies  correspond 
exactly  to  the  peaks  shown  in  Figure  3.9.  The  equivalent  system  is  shown  in  Figure  3.10. 
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Frequency  (Hz) 

Figure  3.9  H  (Q)  for  Free-Free  Beam 
(a-setDOFs:  1, 17) 


Figure  3.10  ABC  Configuration  -  a-set:  [1  17] 
(Measured  coordinates  restrained  to  ground) 
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IV.  ABC  CONFIGURATION  IN  SEN STI VITY -BASED  UPDATING 


A.  SENSITIVITY  MATRIX  DEFINED 

As  previously  mentioned,  the  disparity  in  the  number  of  known  parameters 
(measured  modal  parameters)  versus  the  number  of  parameters  to  be  adjusted  in  order  to 
update  a  FE  model  defines  an  underdetermined  problem.  Another  significant  difficulty  is 
that  the  parameters  that  are  in  error  are  unknown.  Sensitivity-based  updating  is  used  for 
error  localization  in  order  to  find  those  parameters  that  require  adjustment.  The  ABC 
configuration  frequencies  can  be  used  in  addition  to,  or  instead  of,  the  standard  baseline 
system  mode  frequencies  in  model  updating  and  damage  detection.  The  ABC 
frequencies  correspond  to  the  same  structural  system  as  the  baseline  frequencies,  but  with 
different  boundary  conditions.  The  governing  equation  for  sensitivity-based  updating  is 
{A©}  =  [T]{ADV}  (4.1) 

where  {A©}  is  a  vector  of  natural  frequency  errors.  The  errors  are  the  difference  between 
the  experimental  natural  frequencies  and  the  analytical  natural  frequencies  {©x-  ©a}. 

The  {ADV}  term  is  the  vector  of  changes  to  be  calculated  for  specified  model 
parameters,  known  as  "design  variables,"  and  [T]  is  the  sensitivity  matrix.  This  matrix  is 
solved  using  the  eigensensitivity  method  shown  in  Ref.  [6].  [T]  is  defined  as: 
{<|>baseiine}T[Ak]  { baseline},  where  Ak  =  [kx  -  ka]  and  [T]  is  of  size  m  x  p  where  m  is  the 
number  of  modes  and  p  is  the  number  of  design  variables  used.  Each  ABC  system 
defines  additional  rows  of  Eq.  (4.1),  defining  the  equation, 

{A©'}  =  [T']{ADV}  (4.2) 
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where  Tk  =  dcof  /dDVJ ,  and  cok  is  the  i*  natural  frequency  of  the  kth  ABC  configuration 

system.  Baseline  system  quantities  are  identified  with  the  superscript  "0"  and  ABC 
systems  are  identified  with  superscripts  from  1  to  "k"  where  k  is  the  total  number  of  ABC 


systems  used.  Combining  the  equations, 


fK)l 

iA“'} 

>  =  < 

[T’l, 

.VI: 

.[1*1. 

{ADV} 


(4.3) 


The  degree  of  coupling  between  the  ABC  systems  and  the  baseline  system  can  be 


adjusted  by  deleting  or  retaining  individual  columns  of  [Tk].  Partial  coupling  can  be 
established  by  partitioning  such  that  some  of  the  DV's  are  associated  with  the  baseline 
system  only,  some  with  both  the  baseline  and  the  ABC  system,  and  some  with  the  ABC 
system  only.  [Ref.  1] 

Using  the  ABC  system  sensitivities  helps  to  eliminate  or  greatly  reduce  the 
problem  of  poorly  conditioned  or  rank  deficient  [T]  matricies.  Columns  of  [T°]  can  be 
replaced  with  columns  of  [Tk]  in  order  to  improve  conditioning.  Two  closely  spaced 
elements  in  a  model  will  always  have  columns  of  [T°]  which  are  nearly  dependent.  This 
prevents  the  localization  of  the  error  and  discrimination  of  the  damage  between  the 
elements.  By  replacing  a  column  of  [T°]  with  a  column  of  [Tk],  the  columns  are  no 
longer  dependent  and  error  localization  and  damage  detection  is  improved,  resulting  in  a 
more  accurately  updated  FE  model. 
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B.  IMPROVED  CONDITIONING  USING  ABC  SENSITIVITIES 

Applying  the  case  of  a  free-free  beam,  the  sensitivity  matrix  [T]  is  generated  for 
the  undamaged  structure  model  in  order  to  solve  for  the  changes  in  the  element  El  values 
representing  damage,  i.e.  AEI  =  ADV,  and 

{ffl0}  =  [T°]{AEI}  (4-4) 

To  demonstrate,  sensitivities  will  be  calculated  for  all  10  elements  and  for  the  first  10 
elastic  modes  resulting  in  a  [T°]  matrix  of  size  10  by  10.  A  calculation  of  the  rank  of  this 
matrix  reveals  that  it  is,  in  fact,  rank  deficient:  Rank(T°)  =  5.  This  will  not  provide  a 
fully  determined  solution  for  {AEI}.  With  excitation  applied  at  DOFs  1  and  17,  it  is 
possible  to  identify  the  ABC  system  frequencies  with  these  DOF  constrained  to  ground. 
(See  Figures  3.9  and  3.10)  The  corresponding  [T1]  matrix  is  generated  and  is  also  of  size 
10  by  10, 

{©^[T’KAEI}  (4-5) 

with  the  calculated  rank:  Rank(T’)  =  10. 

The  [T1]  matrix  is  full  rank  due  to  the  asymmetric  ABC  configuration. 

C.  COMPUTER  SIMULATED  DAMAGE  DETECTION  USING  ABC 
SYSTEM  SENSITIVITIES 

Four  examples  have  been  constructed  to  demonstrate  the  usefulness  of  a  single 
ABC  system  set  of  frequencies  and  sensitivity  matrix  in  lieu  of  the  standard  system  set. 
The  ABC  system  described  by  Figures  3.9  and  3.10,  a-set:  [1  17],  is  again  used.  In  these 
four  examples,  changes  in  elemental  El’s  have  been  user  inputted  to  simulate  a  damaged 
structure.  Two  subsequent  examples  will  show  the  usefulness  of  this  technique  when 
damage  location  and  magnitude  is  unknown. 
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1. 


Damage  Detection:  Example  1 


For  this  example,  a  10%  reduction  in  El  is  made  at  elements  #2  and  #3.  The 
natural  frequencies  of  the  undamaged  (FE/Analytic)  and  simulated  damaged 
(Test/Experimental)  models  are  shown  in  Table  4.1 .  Using  the  baseline  system,  the 
largest  full  rank  submatrix  of  [T°]  is  5  by  1 0,  and  the  best  solution  for  the  damage  was  the 
solution  using  modes  1  to  4.  This  solution  is  shown  in  Figure  4.1(a).  The  height  of  the 


Mode 

FE/Analytic 

Test 

%  Error 

Number 

(Hz) 

(Hz) _ 

1 

31.30 

31.13 

2 

86.30 

84.99 

1.5 

3 

169.29 

166.37 

1.7 

4 

280.26 

276.61 

1.3 

5 

419.78 

414.52 

1.3 

6 

588.73 

581.00 

1.3 

7 

788.17 

779.81 

1.1 

8 

1018.0 

9 

1264.6 

10 

1679.4 

1662.9 

0.9 

Table  4.1  Example  1:  FE  and  Test  Frequencies  and  Error 


bars  in  the  plot  indicates  the  magnitude  of  the  error  for  each  element.  In  this  case,  we 
know  that  the  exact  solution  is  10%  in  both  element  #2  and  element  #3.  If  the  ABC 
system  sensitivities  [T1]  and  mode  frequencies  are  used,  the  largest  full  rank  submatrix  is 
10  by  10,  and  the  best  solution  for  the  damage  was  the  solution  using  all  10  modes.  This 
solution  is  shown  in  Figure  4.1  (b).  This  was  chosen  as  the  first  example  because  in  this 
case,  the  baseline  system  provided  a  reliable  result.  It  should  be  noted  that  the  ABC 
system  also  provided  reasonable  results  (other  elements  were  affected  by  1%  or  less).  It 
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will  be  seen  in  the  following  examples  that  the  baseline  system  rarely  provides  reliable 
results  and  in  many  circumstances  actually  provides  results  in  misleading  information. 


0.2 
0.15 
0.1 
0.05 

°  1  23456789  10 

(a) 

0.2 
0.15 
0.1 
0.05 

°  123456789  10 

(b) 

Figure  4.1  Damage  Detection  —  Example  1 
(a)  Baseline  (b)  ABC  system 

2.  Damage  Detection:  Example  2 

For  this  example  a  10%  reduction  in  El  is  made  at  element  #4  and  a  15% 
reduction  in  El  is  made  at  element  #6.  The  natural  frequencies  of  the  undamaged 
(FE/ Analytic)  and  simulated  damaged  (Test)  models  are  shown  in  Table  4.2.  Using  the 
baseline  system,  two  solutions  seemed  to  localize  the  error.  The  first 
was  the  solution  using  modes  1  to  5.  This  solution  is  plotted  in  Figure  4.2(a).  The 
second  was  the  solution  using  modes  1  to  4  and  this  solution  is  shown  in  Figure  4.2(b). 
However,  it  was  known  from  the  start  that  the  damage  was  located  at  elements  4  and  6. 
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Mode 

FE/Analytic 

Test 

%  Error 

Number 

(Hz) 

1 

31.30 

2 

86.30 

85.16 

1.3 

3 

169.29 

167.19 

1.2 

4 

280.26 

276.22 

1.4 

5 

419.78 

413.91 

1.4 

6 

588.73 

1.3 

7 

788.17 

1.3 

8 

1018.0 

1004.5 

1.3 

9 

1264.6 

1246.0 

1.5 

10 

1679.4 

1656.7 

1.4 

Table  4.2  Example  2;  FE  and  Test  Frequencies  and  Error 
In  a  standard  test,  there  would  be  no  prior  knowledge  of  damage  and  the  results  using  the 
baseline  system  would  appear  to  be  accurate.  This  misleading  information  would  not 
help  in  the  updating  of  the  FE  model.  Using  the  ABC  system  and  modes  1  through  8,  the 
solution  shown  in  Figure  4.2(c)  is  obtained  and  proves  to  be  a  reliable  solution  based  on 

the  known  error. 

0.2 

£0.1 
0 

0.2 
g  0.1 
0 

0.2 

u  0.1 

1  2  3  4  5  6  7  8  9  10 


Figure  4.2  Damage  Detection  -  Example  2 
(a)  &  (b)  Baseline  (c)  ABC  system 
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3. 


Damage  Detection:  Example  3 


For  Example  3,  two  elements  at  opposite  ends  of  the  beam  were  chosen.  A  15% 
reduction  in  El  was  made  at  element  #1  and  a  20%  reduction  in  El  was  made  at  element 
#10.  The  natural  frequencies  of  the  undamaged  and  simulated  damaged  models  are 
shown  in  Table  4.3.  Using  the  baseline  system,  the  only  solution  that  provided  any 


Mode 

Number 

FE/Analytic 

(Hz) 

Test 

(Hz) 

%  Error 

1 

31.30 

31.29  ! 

0.0 

2 

86.30 

86.21 

0.1 

3 

169.29  ~ 

168.78 

0.3 

4 

280.26 

278.44 

0.6 

5 

419.78 

415.12 

1.1 

6 

588.73 

579.27 

1.6  1 

7 

788.17 

772.08 

^  2.0 

8 

1018.0 

994.76 

2.3 

9 

1264.6 

1240.1 

1.9 

10 

'  1679.4 

1646.6 

2.0  1 

Table  4.3  Example  3:  FE  and  Test  Frequencies  and  Error 


error  localization  was  the  solution  using  modes  1  to  5.  This  solution  is  shown  in  Figure 
4.3(a).  It  is  immediately  obvious  that  this  solution  is  incorrect,  as  it  indicates  an  error  of 
approximately  42%  in  element  #10  and  no  damage  in  any  other  elements.  Using  the 
ABC  system  and  modes  1  through  7,  an  improved  solution  is  obtained  and  is  shown  in 
Figure  4.3(b). 

4.  Damage  Detection:  Example  4 

For  the  fourth  example  three  elements  were  chosen  to  have  reduced  El  values.  A 
10%  reduction  in  El  is  made  at  element  #3,  a  15%  reduction  at  element  #4,  and  a  20% 
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0.5 


0.4 

0.3 

0.2 

0.1 

0 

0.25 

0.2 

0.15 

0.1 

0.05 

0 


reduction  at  element  #5.  The  natural  frequencies  of  the  undamaged  and  simulated 
damaged  models  are  shown  in  Table  4.4.  The  best  baseline  system  solution  was  obtained 
using  modes  1  to  3  and  is  displayed  in  Figure  4.4(a).  The  solution  localized  error  in  three 
elements,  but  localized  damage  in  element  #7  instead  of  element  #4.  This  is  obviously 
incorrect.  Using  the  ABC  system  and  the  same  modes  (1  to  3),  the  solution  displayed  in 
Figure  4.4(b)  was  obtained.  This  solution  also  is  incorrect  and  misleading.  By  using  the 
ABC  system  and  modes  1  through  6,  however,  an  improved  and  correct  solution  is 
obtained  and  is  shown  in  Figure  4.4(c).  It  is  important  to  note  that  in  this  fourth  example, 
the  ABC  system  also  provided  one  misleading  solution.  Given  the  number  of  modes 


1  23456789  10 


(a) 


1  23456789  10 


(b) 

Figure  4.3  Damage  Detection  -  Example  3 
(a)  Baseline  (b)  ABC  system 
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Mode 

FE/Analytic 

Test 

%  Error 

Number 

(Hz) 

(Hz) 

1 

31.30 

29.87 

4.6 

2 

86.30 

83.84 

2.9 

3 

169.29 

164.68 

2.7 

4 

280.26 

273.90 

2.3 

5 

419.78 

409.57 

2.4  1 

6 

588.73 

574.71 

2.4 

7 

788.17 

768.17 

2.5 

8 

1018.0 

993.73 

2.4 

9 

1264.6 

1231.7 

2.6  i 

10 

1679.4 

1637.8 

2.5 

Table  4.4  Example  4:  FE  and  Test  Frequencies  and  Error 


used,  the  solution  using  the  first  six  modes  would  be  considered  more  reliable  even  if  the 
damage  had  not  been  known  a  priori  due  to  its  improved  conditioning. 


Figure  4.4  Damage  Detection  -  Example  4 
(a)  &  (b)  Baseline  (c)  ABC  system 
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For  the  final  two  examples,  a  random  number  generator  is  used  to  place  errors  in 
1  to  2  elements  of  the  model  by  reducing  El  in  each  of  those  elements  by  up  to  50%. 

Using  this  method,  the  damage  location  and  magnitude  is  unknown  and  left  to  be 
determined  by  using  both  baseline  and  ABC  system  sensitivities.  Examining  the  element 
El  matrix  from  the  MATLAB  code  checks  the  results  determined  by  the  sensitivities. 

5.  Damage  Detection:  Example  5 

The  damage  detection  results  obtained  by  using  the  maximum  number  of  modes 
possible  in  the  baseline  system ,  five  (rank=5),  are  shown  in  Figure  4.5  (a).  It  appears 
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from  this  plot  that  element  #6  is  the  damaged  element.  However,  when  modes  1  through 
4  are  used  in  the  baseline  system  (Figure  4.5(b))  it  appears  that  #5  is  the  “damaged” 
element.  By  using  modes  1  through  5  in  the  ABC  system,  Figure  4.5(c)  was  obtained  and 
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Figure  4.5  Damage  Detection  -  Example  5 
(a)  &  (b)  Baseline  (c)  ABC  system 
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verified  that  #5  was  the  damaged  element.  This  was  believed  to  be  the  best  localization  of 
error  and  the  element  El  matrix  was  checked  for  verification.  The  only  damaged  element 
was  element  #5  and  the  difference  between  damaged  and  undamaged  was  a  23  percent 
reduction,  which  corresponds  to  the  values  in  Figure  4.5(b)  and  (c).  Once  again,  the  use 
of  the  ABC  system  provided  reliable  results  when,  without  the  ABC  configuration,  the 
baseline  information  alone  may  have  been  misleading. 

6.  Damage  Detection:  Example  6 

Once  again,  the  results  using  the  first  five  modes  with  baseline  sensitivity  are 
displayed  in  Figure  4.6(a).  Using  this  solution,  it  appears  that  the  damage  is  located  in 
elements  #2  and  #6.  With  the  first  four  modes,  the  solution  shows  damage  in  elements  #2 
and  #5  as  shown  in  Figure  4.6(b).  Figure  4.6(c)  shows  the  results  obtained  when  using 
all  10  modes  and  the  ABC  system  sensitivity.  It  confirms  that  elements  #2  and  #5  are  the 
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Figure  4.6  Damage  Detection  -  Example  6 
(a)  &  (b)  Baseline  (c)  ABC  system 
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damaged  elements.  Checking  the  El  element  matrix,  it  was  confirmed  that  element  #2 
had  a  10%  reduction  in  El  and  element  #5  had  a  3%  reduction  in  EI.  These  values 
correspond  to  those  seen  in  Figure  4.6(b)  and  (c).  In  summary,  all  6  examples  have 
provided  valuable  information  concerning  damage  detection.  Many  times,  a  reasonable 
solution  can  be  obtained  using  baseline  sensitivities.  However,  with  this  always  seems  to 
be  some  misleading  information  when  the  baseline  solutions  are  used.  The  ABC  system 
will  most  always  provide  the  proper  solution,  but  is  not  perfect  either.  When  the  ABC 
system  is  used  in  conjunction  with  the  baseline  system,  improved  and  accurate  solutions 
are  obtained.  Nowhere  is  this  clearer  than  in  the  previous  two  examples. 
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V.  EXPERIMENTAL  APPLICATION 


So  far,  the  use  of  the  ABC  configuration  to  determine  natural  frequencies  of  a 
system  under  a  variety  of  boundary  conditions  has  been  demonstrated  using  MATL  AB 
and  computer  generated  finite  elements.  To  fully  examine  this  process,  a  laboratory 
experiment  was  constructed. 

A.  EQUIPMENT  SETUP 

A  mild  steel  beam,  60.75  inches  in  length,  1.5  inches  wide  and  0.5  inches  thick, 
was  suspended  from  the  overhead  by  monofilament  fishing  line  in  order  to  represent  a 
free-free  beam  with  two  rigid  body  modes.  The  beam  was  instrumented  with  5  PCB 
Piezotronics  Inc.  accelerometers,  each  15.25  inches  apart.  Each  of  the  accelerometers 
was  wired  into  a  5  kHz  Zonic  System  7000  front  end  digital  signal  processor  (DSP) 
consisting  of  2  digital  chassis  with  16  channels  each  and  1  HP  workstation.  An  excitation 
was  applied  successively  to  the  opposite  side  of  each  of  the  accelerometers  by  a  25  lb 
peak  force  shaker  and  lateral  excitation  stand  Model  2050A,  manufactured  by  The  Modal 
Shop,  Inc.  A  MB  Dynamics  Model  SS250VCF  amplifier  set  to  constant  voltage  powered 
the  shaker.  The  shaker  was  operated  at  0.25  volts  (rms).  The  applied  force  was  wired 
into  the  DSP  as  Channel  1,  and  the  5  accelerometers  as  Channels  2  through  6, 
respectively.  The  Zonic  system  was  connected  to  the  I-DEAS  Master  Series  5  software 
operated  in  the  Signal  Processing  mode.  The  laboratory  setup  is  shown  in  Figure  5.1. 
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Figure  5.1  Laboratory  Setup 


B.  DATA  COLLECTION 


Using  I-DEAS,  FRFs  were  collected  when  the  force  was  applied  at  each  node. 

For  each  reference  (i.e.  for  the  force  applied  at  each  node),  the  response  at  each  node  was 
measured.  I-DEAS  was  configured  to  measure  1601  spectral  lines  over  frequencies 
ranging  from  0  -2000  Hz.  Each  shaker  provided  a  1601  by  5  matrix  containing  the 
complex  amplitudes  every  1.25  Hz  for  a  given  excitation.  The  combined  result  provided  a 
1601  by  25  matrix  for  the  entire  system.  This  matrix  was  exported  to  MATLAB  in  order 
to  compute  the  o-set  frequencies  of  the  system. 

C.  O-SET  FREQUENCIES  USING  EXPERIMENTAL  DATA 

In  order  to  establish  a  comparison  and  validate  the  results  obtained  in  the 
laboratory,  a  finite  element  model  of  the  beam  was  created.  The  FE  model  contained  the 
given  dimensions  of  the  beam  and  a  Modulus  of  Elasticity  (E)  of  27.3  x  10  lbf/in  , 
corresponding  to  mild  steel.  This  E  value  was  chosen  in  order  to  bring  the  baseline 
model  frequencies  closer  to  the  test  frequencies.  It  was  composed  of  64  elements  with  65 
nodes  (2  DOF  per  node  for  130  DOF  total).  Nodes  1,  17,  33,  49,  and  65  corresponded  to 
transducer  locations  1  through  5,  respectively,  on  the  real  beam. 

The  theories  presented  earlier  and  demonstrated  through  computer  simulation 
were  applied  to  the  experimental  data.  By  using  ABC  configurations,  the  frequencies  for 
the  beam  under  30  different  boundary  conditions  can  be  found  from  the  baseline  (free- 
free)  information  without  having  to  reconfigure  the  experimental  setup.  Table  5.1 
contains  the  first  ten  natural  frequencies  (in  Hz)  of  the  system  computed  from  the  test  and 
compares  them  to  the  natural  frequencies  determined  by  the  FE  model.  Figure  5.2 
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contains  a  plot  of  the  1,1  element  of  the  FRF  (H(Q))  matrix  and  displays  the  first  five 
modes. 


FE/ Analytic  (Hz) 

Test  (Hz) 

%  Difference 

29.1 

29.0 

0.3 

80.3 

79.5 

157.4 

155.3 

1.3 

260.2 

257.8 

0.9 

388.7 

388.3 

0.1 

542.8 

534.0 

1.6 

722.7 

1.9 

928.3 

912.9 

1.7 

1159.6 

1104.1 

4.8 

1416.6 

1393.5 

1.6 

Table  5.1  FE  &  Test  frequencies  for  the  free-free  beam 


Figure  5.2  Hn(Q)  for  free-free  beam 
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The  first  o-set  system  computed  was  the  system  with  node  1  restrained  to  ground, 
corresponding  to  the  system  shown  in  Figure  5.3.  The  FRF  matrix  was  inverted  to  find 
the  impedance  matrix  for  the  a-set  ( H^j  (Q))  and  the  peaks  of  the  impedance  matrix 

correspond  directly  to  the  natural  frequencies  of  the  system  with  node  1  restrained  to 
ground.  These  o-set  frequencies  were  compared  to  the  o-set  frequencies  determined  by 
the  FE  model  and  are  shown  in  Table  5.2.  A  plot  of  (Q)  showing  the  first  five 
frequencies  is  seen  in  Figure  5.4. 

•  •  • .  'll 

Figure  5.3  ABC  configuration:  Node  1  restrained  to  ground 


FE/Analytic  (Hz) 

Test  (Hz) 

%  Difference 

20.0 

20.0 

0.0 

65.0 

65.6 

0.9 

135.7 

136.0 

0.7 

232.1 

235.0 

1.2 

354.1 

356.2 

0.6 

501.9 

503.8 

0.4 

675.3 

669.6 

0.8 

874.5 

873.3 

0.1 

1099.3 

1080.7 

1.7 

1349.9 

1354.1 

0.3 

Table  5.2  FE  &  Test  Frequencies  for  Node  1  restrained  to  ground 
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Figure  5.4  (Q)  for  Node  1  restrained  to  ground 

In  both  cases,  the  test  frequencies  correspond  nicely  with  the  frequencies 
predicted  by  the  FE  model.  This  validates  the  theories  and  expectations  described  in  the 
earlier  chapters.  With  only  free-free  data  available,  the  natural  frequencies  of  the  system 
under  any  boundary  configuration  can  be  found.  The  number  of  transducers  (i.e. 
measurement  locations)  limits  the  number  of  systems  for  which  frequencies  are  available 
and  in  the  given  example,  up  to  30  system  frequencies  can  be  found.  Six  more  ABC 
systems  were  computed  and  the  results  for  the  first  ten  modes  are  tabulated  in  Tables  5.3 
-  5.8.  The  system  configurations  and  H~a'  (Q)  plots  are  shown  in  Figures  5.5  -  5.16. 
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Figure  5.5  ABC  Configuration:  Nodes  1  &  2  restrained  to  ground 


Test  (Hz) 

%  Difference 

6.8 

7.2 

5.9 

44.3  1 

44.7 

0.9 

125.7 

125.9 

0.2  1 

228.9 

231.9 

1.3  1 

293.4 

297.0 

1.3  1 

443.1 

443.9 

0.2 

652.7 

648.4 

h  0.7 

868.4 

867.9 

0.1 

990.2 

980.5 

1.0 

1252.2 

1253.1 

0.1 

Table  5.3  FE  &  Test  Frequencies  for  Nodes  1,2  restrained  to  ground 
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Figure  5.6  H  J  (Q)  for  Nodes  1 ,2  restrained  to  ground 
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Figure  5.7  ABC  Configuration:  Nodes  1,4  restrained  to  ground 


FE/ Analytic  (Hz) 

Test  (Hz) 

%  Difference 

49.4 

49.1 

0.4 

108.0 

110.0 

1.8 

219.5 

223.9 

2.0 

349.6 

347.6 

0.6 

427.6 

427.2 

0.1 

606.6 

618.8 

2.0 

848.7 

856.0 

0.9 

1092.5 

1045.2 

4.3 

1228.7 

1229.7 

0.1 

Table  5.4  FE  &  Test  Frequencies  for  Nodes  1,4  restrained  to  ground 
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Figure  5.8  (Q)  for  Nodes  1 ,4  restrained  to  ground 
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Figure  5.9  ABC  Configuration:  Nodes  1 ,2,3  restrained  to  ground 


FE/Analytic  (Hz) 

Test  (Hz) 

%  Difference 

14.5 

15.0 

3.4 

96.5 

96.8 

0.3 

220.8 

225.1 

1.9 

282.2 

283.4 

0.4 

375.0 

381.1 

1.6 

588.3 

537.4 

* 

00 

853.2 

858.5 

0.6 

970.9 

957.7 

1.4 

1133.3 

1126.9 

0.6 

1489.3 

1463.6 

1.7 

Table  5.5  FE  &  Test  Frequencies  for  Nodes  1,2,3  restrained  to  ground 
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Figure  5.10  (Q)  for  Nodes  1,2,3  restrained  to  ground 
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Figure  5.11  ABC  Configuration:  Nodes  1,2, 3, 4  restrained  to  ground 


FE/Analytic  (Hz) 

Test  (Hz) 

%  Difference 

49.3 

48.3 

2.0 

215.3 

220.0 

2.2 

277.4 

281.2 

1.2 

367.5 

369.5 

0.5 

448.2 

448.1 

0.0 

841.6 

848.2 

0.8 

961.6 

956.5 

0.5 

1123.0 

1123.2 

1260.9 

1258.4 

0.2 

1879.2 

1895.6 

0.9 

Table  5.6  FE  &  Test  Frequencies  for  Nodes  1,2, 3, 4  restrained  to  ground 
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Figure  5.12  Haa  (Q)  for  Nodes  1 ,2,3,4  restrained  to  ground 
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Figure  5.13  ABC  Configuration:  Nodes  1,2, 3, 4, 5  restrained  to  ground 


%  Difference 

210.6 

2.4 

239.8 

245.2 

2.3 

321.1 

328.2 

2.2 

415.0 

422.1 

1.7 

822.3 

831.1 

1.1 

892.4 

894.8 

0.3 

1040.7 

1043.7 

0.3 

HIHIHDSEIsEHHHi 

1214.1 

1.1 

1869.1 

1.0 

1954.5 

1952.8 

0.1 

Table  5.7  FE  &  Test  Frequencies  for  Nodes  1,2, 3 ,4,5  restrained  to  ground 


0  50  100  150  200  250  300  350  400  450  500 


Figure  5.14  Hj  (Q)  for  Nodes  1,2, 3, 4, 5  restrained  to  ground 
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Figure  5.15  ABC  Configuration:  Nodes  1,3,5  restrained  to  ground 


FE/ Analytic  (Hz) 

Test  (Hz) 

%  Difference 

51.4 

52.2 

1.6  ~1 

80.3 

82.4 

2.6 

205.6 

210.7 

2.5 

260.2 

266.5 

2.4 

462.5 

471.1 

1.9 

542.8 

554.8 

2.2 

822.3 

831.8 

1.2 

928.3 

940.9 

1.4 

1284.8 

1306.8 

1.7 

141 6.6 

1428.6 

0.8 

Figure  5.8  FE  &  Test  Frequencies  for  Nodes  1,3,5  restrained  to  ground 
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Figure  5.16  Hj  (Q)  for  Nodes  1,2, 3, 4, 5  restrained  to  ground 


D.  DAMAGE  DETECTION  USING  EXPERIMENTAL  DATA 


The  next  logical  step  to  take  is  to  try  to  localize  error  in  a  FE  model  by  using  the 
experimental  data.  In  Chapter  IV,  sensitivity-base  damage  detection  was  used  to  localize 
error  in  a  "damaged"  FE  model  by  using  another  FE  model  with  no  damage.  We  now 
have  a  set  of  data  for  a  real  beam  that  has  no  damage.  The  experimental  data  is  now 
needed  in  order  to  find  the  error  in  a  FE  model  of  the  same  system.  To  perform  these 
calculations,  an  FE  model  was  again  created,  this  time  of  16  elements.  Fewer  elements 
were  used  in  order  to  decrease  computational  time;  the  actual  change  in  frequencies  from 
those  determined  using  64  elements  was  negligible.  A  20%  reduction  in  El  was  made  at 
element  #8  in  the  FE  model.  A  model  of  the  free-free  beam  (baseline)  was  created  as 
well  as  a  model  of  the  ABC  system  with  nodes  1  and  4  restrained  to  ground. 

To  check  the  method,  FE  frequencies  were  substituted  for  the  test  frequencies  in 
both  the  baseline  system  and  the  ABC  system.  The  results  are  shown  in  Figure  5.17. 


(b) 

Figure  5.17  Error  localization  using  FE  frequencies 
(a)  Baseline  (b)  ABC  system 
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The  error  was  localized  using  the  first  two  modes  of  the  baseline  system  (Figure  5. 1 7(a)) 
and  then  was  confirmed  using  the  first  two  modes  of  the  ABC  system  (Figure  5.17(b)). 
This  validated  the  method  and  now  the  FE  frequencies  were  replaced  by  the  actual  test 
frequencies  for  the  system.  These  results  are  shown  in  Figure  5.18. 


Figure  5.18  Error  localization  using  test  frequencies 
(a)  Baseline  (b)  ABC  system 


The  first  two  modes  were  again  used  for  both  the  baseline  system  (a)  and  the 
ABC  system  (b)  since  these  were  the  modes  that  resolved  the  error  when  the  FE 
frequencies  were  used.  This  time,  however,  the  baseline  system  localized  a  20%  error  in 
element  #8  as  well  as  a  13%  error  in  element  #12.  Clearly  this  is  incorrect.  The  ABC 
system  provided  better  results,  but  also  indicated  error  in  element  #12  at  5%.  The 
differences  in  the  frequencies  between  the  test  data  and  the  FE  data,  however  slight,  are 
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the  sources  of  the  discrepancies.  In  the  first  two  modes  of  the  free-free  beam,  the  test 
frequencies  differed  from  the  FE  frequencies  by  0.3%  and  1.0%,  respectively  (see  Table 
5.1).  For  the  ABC  system’s  first  two  modes,  the  test  frequencies  differed  by  0.0%  and 
0.4%,  respectively  (see  Table  5.4).  Therefore,  in  order  to  detect  and  localize  errors  in  FE 
models  with  test  information,  it  is  imperative  that  the  error  between  test  and  FE 
frequencies  be  minimized.  Suggested  approaches  on  how  to  achieve  this  are  proposed  in 
the  Recommendations. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  main  objective  of  this  thesis  was  to  examine  a  method  in  which  the  natural 
frequencies  of  a  system  under  a  variety  of  boundary  conditions  could  be  found  by  using 
only  one  experimental  database.  Previously,  experiments  had  to  be  reconfigured  in  order 
to  find  more  than  the  baseline  frequencies.  This  results  in  more  time  and  inevitably,  an 
increase  in  cost  for  structural  analysis.  Secondly,  a  method  by  which  model  error 
detection  and  localization  could  be  improved  was  investigated. 

A.  CONCLUSIONS 

The  following  conclusions  can  be  drawn  form  the  analyses  presented: 

1 .  The  natural  frequencies  of  a  structure  under  various  boundary  conditions  are 
available  from  a  baseline  FRF  matrix.  Inverting  a  spatially  incomplete  FRF  at 
each  frequency  results  in  an  impedance  spectrum  whose  peaks  correspond  to 
the  natural  frequencies  of  the  structure  restrained  to  ground  at  measure 
coordinates. 

2.  With  both  the  baseline  and  ABC  frequencies  available,  structural  sensitivities 
can  be  generated  from  both  the  baseline  and  from  the  system  with  boundary 
conditions  artificially  applied  corresponding  to  the  FRF  matrix  (a-set).  The 
addition  of  the  ABC  data  has  proven  to  be  valuable  in  conjunction  with  the 
baseline  information  in  localizing  error  and/or  damage. 
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B.  RECOMMENDATIONS 


The  following  are  areas  of  recommended  improvement  to  continue  with  this  field 
of  study.  The  first  two  are  recommendations  for  improving  sensitivity-based  error 
localization  when  using  test  data. 

1 .  Curve  fit  the  FRF  in  I-DEAS  prior  to  exporting  the  files  into  MATLAB. 

2.  Once  in  MATLAB,  curve  fit  the  impedance  spectra  near  the  ABC  frequency 
to  identify  the  ABC  natural  frequencies. 

3.  Continue  with  the  use  of  ABC  systems  on  more  advanced  structures  to  further 
this  study. 
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APPENDIX  MATLABCODE 


A.  MAIN  PROGRAMS 


Build2Beams.m 


%  This  program  assembles  the  mass  and  stiffness  matrices  for  2  beams,  with  user 
%  defined  boundary  conditions  referred  to  as  "BeamA"  (analysis)  and  "BeamX" 

%  (experimental).  The  program  can  be  run  in  both  “build”  and  “read”  modes.  In  the 
%  “build  mode,”  the  user  provides  baseline  data  for  BeamA,  assumed  to  be  a 

%  homogeneous,  uniform  beam.  Data  provided: 

% 

%  (1)  Beam  length 

%  (2)  Number  of  elements 

%  (3)  Nominal  El 

%  (4)  Nominal  cross-sectional  area 

%  (5)  Nominal  weight  density 

% 

%  The  program  then  prompts  the  user  for  instructions  on  how  to  modify  "BeamA" 

%  data  to  arrive  at  "BeamX"  data.  The  user  can  modify  element  El  values  and  the 

%  modification  can  be  applied  to  either  a  single  element,  or  range  of  elements.  The 

%  beam  definition  data  is  saved  in  a  binary  (.mat)  file  "beamdata"  at  the  end  of 
%  execution. 

% 

%  In  the  "read"  mode,  the  program  will  load  the  beam  data  written  at  the  last 
%  execution  of  the  program  and  will  then  allow  the  user  to  modify  the  data  for 

%  BeamX.  By  answering  "y"  to  this  prompt,  the  program  will  copy  BeamA  data  into 

%  BeamX  data,  and  then  prompt  the  user  (identically  as  described  above)  for  data  on 

%  how  to  adjust  the  BeamA  data  in  order  to  arrive  at  the  BeamX  data. 

% 

%  The  program  calls  the  following  support  programs: 

% 

%  BeamA_Prompt.m 

% 

%  BeamX_Prompt.m 

%  BoundaryConditions.m 
%  Assemble2Beams.m 

% 

% 

%  clearspace.m 

% 

%  Written  by  J.H.  Gordis 

%  Modified  by  C.P.  DeGregory 


--  Prompts  User  for  BeamA  nominal  beam  data 
(“build”  mode  only) 

--  Prompts  User  for  BeamX  modification  beam  data 

—  Prompts  user  for  B.C.'s  and  applies  them. 

—  Called  by  Build2Beams,  builds  [ka]  [ma]  [kx]  [mx], 
plots  freqs. 

Saves  data  to  "beamdata.mat" 

—  Clears  all  unnecessary  variables  from  workspace. 
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disp( ' '); 

disp( '  Build2Beams  execution...'); 


input_flag  =  'b'; 

input_flag  =  input('  "B"uild  new  beams  or  "R"ead  from  file  ?  (b/r)  ','s'); 

%  Build_or_read  if  loop: 

0^  ******************************************************************** 

if  input_flag  ==  'b';  %  (IF#1 )  Prompt  user  for  beam  definition  data 

0^  ******************************************************************** 


dispC  Building  2  beams  from  scratch...') 

BeamA_Prompt;  %  Prompt  for  BeamA  Data:  run  prompt  script 

BeamX_Prompt;  %  Prompt  for  BeamX  Modification  Data:  run  prompt  script 

Assemble2Beams;  %  Run  script  to  assemble  mass  and  stiffness. 
BoundaryConditions;  %  Prompt  for,  and  apply  boundary  conditions 

%  ***  Save  Defining  Parameters  for  Beams  *** 

disp('  ...saving  beam  data  to  file') 

save  beamdata  connect  ndof  element_length  element_EI  element_mass  mass_diag 
restraint_switch  free_dof_set... 
num_elements  numrbm  lama  lamx 
%save  beamdata  phia  phix 

0/Q  ********************************************************************* 
elseif  input_flag  =  'r';  %  (IF#1)  file  read  for  beam  definition  data 

%************ *********************************************** *********** 

disp('  Loading  2  beams  from  ascii  file  "beamdata"...') 
load  beamdata 

%  Prompt  user  to  either  build  two  beams  from  data  in  file,  or  modify  beamx  data: 
modify_beamx  =  'n'; 

modify_beamx  =  input('  Modify  beamX  data  (y/n)?  ','s'); 

if  modify_beamx  ==  'y';  %  (IF#2) 
disp('  Modify  BeamX...') 
reset_beamx  =  'y'; 

reset_beamx  =  input('  Reset  BeamX  properties  to  BeamA  properties?  . . . 
(y/n)  ','s'); 

if  reset_beamx  ==  'y';  %  Copy  beamA  data  to  beamX  data: 
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element_EI(:,2)=element_EI(:,  1 ); 
element_mass(:,2)=element_mass(:,l); 

%  ...Display  BeamA  properties  to  user: 

sprintf('Number  of  elements:  %3i',num_elements) 
end;  %  end  reset_beamx 

BeamX_Prompt;  %  Prompt  for  BeamX  Modification  Data:  run  script 

%  End  modify_beam  if  block  (IF#2) 
end 

Assemble2Beams;  %  Run  script  to  assemble  mass  and  stiffness. 
BoundaryConditions;  %  Apply  boundary  conditions  (that  already  exist) 

%  ***  Save  Defining  Parameters  for  Beams  *** 

disp('  ...saving  beam  data  to  file') 

save  beamdata  num_elements  connect  ndof  element_length  element_EI ... 

element_mass  mass_diag  restraint_switch  ffee_dof_set... 
lama  phia  lamx  phix  num_elements  connect  ka  ma  kx  mx  num  rbm 

0^  ********************************************************************* 
end;  %  End  Build  or  read  if  loop  IF#  1 

0^*  5|e  *******************************************************************  * 

disp('  Build2Beams  end.') 

% 
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BeamSensitivity.m 


%  This  program  calculates  the  stiffness  and  mass  matrix  partials  by  finite 
%  differences.  That  is,  for  example,  the  [k]  matrix  is  assembled  twice,  once  for  the 
%  nominal  beam  parameters,  and  a  second  matrix  is  assembled  based  on  a  small 
%  perturbation  of  element  EL 

% 

%  This  program  makes  use  of  the  beam  data  created  by  the  program 

%  “Build2Beams.m,”  reads  the  data  file  BeamData  and  resets  BeamX  El  data  to  be 

%  the  same  as  BeamA  data. 

% 

%  The  user  is  prompted  for  element  label  numbers  to  be  used  as  design  variables  by 
%  the  script  “BeamSens_Prompt”  which  stores  this  element  label  data  in  the  matrix 
%  “EI_el_lbl.” 

% 

%  EI_el_lbl:  This  array  is  (max_num_dv,num_elements),  where  max_num_dv  is 
%  the  maximum  number  of  design  variables  allowed,  and  each  design  variable  may 
%  include  from  1  to  "num_elements"  elements.  Therefore,  the  sets  of  elements 
%  included  in  two  design  variables  may  share  common  elements. 

% 

%  Based  on  the  number  of  non-zero  rows  in  the  array  “EI_el_lbl,”  matrix  partials 
%  are  calculated  for  each  non-zero  row.  For  example,  to  start,  the  first  row  of 

%  El  el  lbl  is  read,  and  a  small  perturbation  in  element  El  is  applied  simultaneously 

%  to  the  elements  listed  in  this  row. 

% 

%  The  script  generates  a  matrix  containing  mode  frequency  sensitivities: 

% 

%  [sens_EI]:  num_modes  *  size(EI_el_lbl,l) 

% 

% 

%  Written  by  J.H.  Gordis 

o/0  Modified  by  C.P.  DeGregory 

% 

%********************************************************************** 
%  Start  code: 


format  long; 
disp(' '); 

disp('  BeamSensitivity  execution...');disp(' '); 

disp('  Loading  2  beams  from  ascii  file  "beamdata"...');disp(' '); 

load  beamdata 

num  elements  =  num  elements 
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clear  num_mass_dv  num_EI_dv  defme_M_K_dv  EI_el_lbl  mass_el_lbl 
clear  sens_EI  sens  mass 

BeamSens_Prompt;  %  Prompt  for  DV  element  label  definitions 


o/o  *************************  initialization  *************************** 
0,/  ******************************************************************** 


mass_change  =  1 ;  %  Percentage  mass  change  for  sensitivity  calculation. 

EI_change  =  1 ;  %  Percentage  El  change  for  sensitivity  calculation. 

element_mass_orig  =  element_mass;  %  Copy  properties  over  to  retain  them. 
element_EI_orig  =  element_EI;  % 

icnt_dv  =0;  %  Initialize  number  of  DV's  to  zero.  Will  count  up  all  DV's 

%  Prompt  for  number  of  mode  frequencies  to  generate  sensitivities  for: 

% - - - - - 

num_modes  =  input('  Enter  number  of  elastic  modes  for  sensitivity  calculations»  '); 
start_mode  =  num_rbm  +  1 ;  %  Skip  the  rigid  body  modes. 


0,/  ******************************************************************** 
%  **************  EI  SENSITIVITY  CALCULATION  LOOP  ****************** 

o/0  ******************************************************************** 

num_EI_dv  =  size(EI_el_lbl,l); 

for  icnt_dv  =  1  :num_EI_dv;  %  Repeat  all  of  this  for  all  EI  DV's  defined. 

%  Resetting  BeamX  properties  to  BeamA  properties.. .');disp(' '); 
element_EI(:,2)  =  element_EI(:,l); 

%  Get  element  labels  for  current  DV : 

Ibis  =  EI_el_lbl(icnt_dv,find(EI_el_lbl(icnt_dv,:))); 

%  Perturb  EI  props  for  all  elements  in  "Ibis": 

element_EI(lbls,2)  =  element_EI(lbls,2)  *  (1  +  (EI_change/100) ); 

Assemble2Beams;  %  Run  script  to  assemble  baseline  and  perturbed  beams. 
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BoundaryConditions;  %  Apply  boundary  conditions. 

[lam_base,phi_base]=fmodes2(ka,ma);  %  Get  modes  of  baseline  beam 

%  Form  El  derivative  matrices: 

k_delta  =  (kx  -  ka)/(EI_change/l  00);  %  in  %/l  00 

%  Mode  ffeq  sens  loop: 

end_mode  =  start_mode  +  (num_modes  - 1);  %  Get  requested  number  of  modes. 
row_num  =  0; 

for  icnt_modes  =  start_mode:end_mode; 
rowjnum  =  rownum  +1 ; 

sens_EI(row_num,icnt_dv)  =  phi_base(:,icnt_modes)'  *... 
k_delta  *  phi_base(:,icnt_modes); 

end; 


end; 

%  End  "for  icnt_dv"  outer  loop  for  sensitivity  calculations 
%  Run  script  to  clear  out  variables: 

clearspace 

%  Copy  element  El  and  mass  properties  back  into  arrays: 
element_EI  =  element_EI_orig; 
element_mass  =  element_mass_orig; 

save  beamdata  element_EI  element_mass  sens_EI 

%  Save  all  variables  to  beamdata  (essentially  does  an  append  to  existing  beamdata) 


disp(’  BeamSensitivity  end.') 

% _ 
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Beamex.m 


%  This  program  calculates  an  FRF  matrix  and  the  inverse  of  an  FRF  matrix  in  order 

%  to  identify  the  Artificial  Boundary  Condition  frequencies.  It  loads  the  binary 

%  (.mat)  file  “Beamdata”  and  determines  the  natural  frequencies  of  the  system  from 
%  the  stiffness  and  mass  matricies  by  calling  the  function  “fmodes.”  The  program 

%  also  calls  the  function  “fOset_from_Aset”  to  compute  the  o-set  from  a  user 

%  defined  a-set,  as  well  as  “fmodes”  again,  this  time  to  solve  for  the  o-set  natural 
%  frequencies  (which  correspond  to  the  ABC  frequencies).  The  FRF  is 

%  calculated  as  the  inverse  of  the  Impedance  (H  =  Z-l,  where  Z  =  k  -  (2m).  Both 

%  the  FRF  and  inverse  FRF  (Impedance)  plots  are  displayed  over  the  frequency 
%  range  0  -  800  Hz. 


dispC - ') 

load  beamdata 

[lam, phi]  =  fmodes(kx,mx); 
freq  =  sqrt(lam)/2/pi; 

%  User  is  asked  to  input  a-set 
aset  =  input('Enter  a-set  (vector  form): '); 

[oset]  =  fOset_from_Aset(length(mx),aset); 

ko  =  kx(oset,oset); 
mo  =  mx(oset,oset); 

disp(' A-Set  Defined  as:') 
aset 

disp('0-Set  Frequencies') 

[lamo,phio]  =  fmodes(ko.mo); 
freqo  =  sqrt(lamo)/2/pi; 
icnt  =  0; 

freq  =[1:1:  5030]; 
ffeqp  =  ifreq  /  2  /pi; 

for  ifreq  =  ifreq; 


icnt  =  icnt +1; 

h  =  inv(kx  +  .001*sqrt(-l)  *  kx  -  ifreqA2  *  mx); 
hx  =  h(aset,aset); 


55 


hxp  1 1  (icnt)  =  hx(  1 , 1 ); 

%  hxp22(icnt)  =  hx(2,2); 

%  hxp33(icnt)  =  hx(3,3); 

%  hxp44(icnt)  =  hx(4,4); 

%  hxp55(icnt)  =  hx(5,5); 

%  hxp66(icnt)  =  hx(6,6); 

zabc  =  inv(hx); 
zabcl  l(icnt)  =  zabc(l,l); 

%  zabc22(icnt)  =  zabc(2,2); 

%  zabc33(icnt)  =  zabc(3,3); 

%  zabc44(icnt)  =  zabc(4,4); 

%  zabc55(icnt)  =  zabc(5,5); 

%  zabc66(icnt)  =  zabc  (6,6); 


%  Plot  FRF  and  Impedance 
elf 

figure(l) 

plot(freqp,loglO(hxpl  l));grid 
xlabel('Frequency  (Hz)') 
axis([0,800,-7,3]) 

figure(gcf+l) 

plot(freqp,loglO(zabcl  l));grid 
xlabel('Frequency  (Hz)') 
axis([0,800,-l,7]) 

save  hdata  freqp  zabcl  1 
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Solveos.m 


%  This  program  solves  for  the  vector  {(El}  following  user  guidelines.  It  first  loads 
%  data  which  was  saved  after  two  runs  of  the  “Build2Beams”  and 
%  “Beam_Sensitivity”  programs.  The  first  run  is  for  a  baseline  free-free  beam  with 
%  given  element  El  modifications.  The  analytical  and  experimental  (simulated 
%  damaged)  eigenvalues  (lama  &  lamx)  and  the  sensitivity  matrix  (sens_EI)  are 

%  stored  in  beamO.mat.  The  second  run  is  for  a  beam  with  user  defined  constraints 

%  at  given  measurement  locations  (a-set)  and  lama,  lamx  and  sens_EI  for  this  case 
%  are  stored  in  beaml  .mat. 

% 

%  Written  by  J.H.  Gordis 

%  Modified  by  C.P.  DeGregory 


load  beamO 

lamaO=lama; 

lamxO=lamx; 

TO  =  sens_EI; 

clear  lama  lamx  sens_EI 

load  beaml 

lamal=lama; 

lamxl=lamx; 

T1  =  sens_EI; 

clear  lama  lamx  sens_EI 

rerun  =  'y'; 
while  rerun  —  'y' 

%  Prompt  for  number  of  Beam  mode  frequencies: 

%  Note:  for  free-free  structures,  only  elastic  modes  are  in  T  matrix 

num_modesO  =  size(T0,l); 
numjmodesl  =  size(Tl,l); 

sprintf('Total  BeamO  mode  sens,  data  avail:  %3i',  num_modesO) 
sprintf('Total  Beaml  mode  sens,  data  avail:  %3i',  numjmodesl) 

disp(' ');  disp('  For  BeamO: '); 

use_all_modes  =  input('  Use  all  BeamO  modes  ?  (y/n) » ’,'s'); 


if  use_all_modes  —  'n'; 
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disp('  Enter  (elastic)  mode  numbers  to  include...') 
dispC  Use  MATLAB  vector  format>  1  3  5:7  9  ') 
use_modesO  =  input('  »  ',’s'); 
use_modesO  =  eval(['[’,use_modesO, ']’]); 

else 

use_modesO  =  1  :num_modesO; 

end; 

disp(' ');  disp('  For  Beaml : '); 

use_all_modes  =  input('  Use  all  Beaml  modes  ?  (y/n) »  ','s'); 
if  usealljmodes  —  'n'; 

disp('  Enter  (elastic)  mode  numbers  to  include...') 
disp('  Use  MATLAB  vector  format>  1  3  5:7  9  ') 
use_modesl  =  input('  »  ','s'); 
use_modesl  =  eval(['[',use_modesl, ']']); 

else 

use_modesl  =  l:num_modesl; 

end; 

%  Prompt  for  number  of  DV's: 


num_DV0  =  size(T0,2); 
num_DVl  =  size(Tl,2); 

sprintf('Total  BeamO  DVs  avail:  %3i’,  numJDVO) 
sprintf('Total  Beaml  DVs  avail:  %3i',  num_DVl) 


%  Prompt  user  for  number  of  DVs  to  include: 


disp(' ');  disp('  For  BeamO: '); 

sprintf('Number  of  BeamO  El  DVs:  %3i',num_DV0) 
use_all_dv  =  input('  Use  all  BeamO  DVs  ?  (y/n) »  ','s'); 
if  use_all_dv  ~=  'y'; 

disp('  Enter  DV  label(s)  to  include...') 

disp('  Use  MATLAB  vector  format>  1  3  5:7  9  ') 

usedvO  =  input('  »  ','s'); 

,  use_dvO  =  eval(['[',use_dvO, ']']); 
else 

usedvO  =  1  :num_DV0; 


end 
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disp(' ');  disp('  For  Beaml : '); 
sprintfCNumber  of  Beaml  El  DVs:  %3i',num_DVl) 
use_all_dv  =  input('  Use  all  Beaml  DVs  ?  (y/n) » 
if  use_all_dv  ~=  'y'; 

disp('  Enter  DV  label(s)  to  include...') 
disp('  Use  MATLAB  vector  format>  1  3  5:7  9  ') 
use_dvl  =  input('  »  ','s'); 
use_dvl  =  eval(['[',use_dvl ,']']); 
else 

use_dvl  =  l:num_DVl; 
end 


%  Assemble  Composite  Sensitivity  Matrix: 

format  short 

Tboth  =  [TO(use_modesO,use_dvO);Tl(use_modesl,use_dvl)]; 

%  Calculate  Frequency  Error  Vector  for  the  modes  retained: 

del_ffeq  =  [Iamx0(usejmodes0)-lama0(use_modes0);lamxl(use_modesl> 
lamal  (use_modes  1 )] ; 

%  Solve  system  for  DV  changes: 

sprintf('Min  Dim  of  T:  %3i',min(size(Tboth))) 
sprintf('Rank  of  T:  %3i',rank(Tboth)) 
sprintf('Condition  of  T:  %5.3e',cond(Tboth)) 

delta_dv  =  Tboth  \  del_ffeq 

bar(l :  1 0,abs(delta_dv)),grid 

rerun  =  input( '  Rerun  for  different  case  ?  (y/n) »  ','s'); 
end;  %  end  rerun  while  loop 
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Main.m 


%  This  program  is  similar  to  the  “beamex”  program  in  that  it  computes  FRF  and 

%  Impedance  matrices.  This  script,  however,  is  for  a  spring-mass  system.  In 

%  particular,  this  program  solves  for  the  system  shown  in  Figure  3.1 .  The  stiffness 
%  and  mass  matrices  are  computed  using  the  function  “fspringmass”  and  the  values 

%  for  connectivity,  mass  specification  and  fixity.  The  values  required  for  the  system 

%  shown,  in  addition  to  the  a-set  (DOF  #  1 ),  are  contained  in  the  script  “testdata.” 

%  The  1,1  elements  of  the  FRF  and  Impedance  matrices  are  plotted  versus 

%  frequency. 


clear 

testdata 

[k,m]  =  fspringmass(conn,mspec,fix); 
[lam, phi]  =  fmodes(k,m); 

[oset]  =  fOset_from_Aset(length(m),aset); 

ko  =  k(oset,oset); 
mo  =  m(oset,oset); 

icnt  =  0; 

ifreq  =  [0:.01:4]; 
freqp  =  ifreq/2/pi; 
h  =  zeros(size(k)); 
for  ifreq  =  ifreq 

icnt  =  icnt  +  1 ; 
h  =  inv(k  -  ifreqA2*m); 


hxpll(icnt)  =  h(l,l); 
%hxp22(icnt)  =  h(2,2); 
%hxp33(icnt)  =  h(3,3); 
%hxp44(icnt)  =  h(4,4); 
%hxp55(icnt)  =  h(5,5); 
%hxp66(icnt)  =  h(6,6); 


zabcl  l(icnt)  =  inv(h(l,l)); 
%zabc22(icnt)  =  inv(h(2,2)); 
%zabc33(icnt)  =  inv(h(3,3)); 
%zabc44(icnt)  =  inv(h(4,4)); 
%zabc55(icnt)  =  inv(h(5,5)); 
%zabc66(icnt)  =  inv(h(6,6)); 


end 


60 


orient  landscape 
figure(l) 

plot(freq*2*pi,logl0(hxpl  1)),  grid 
xlabel(‘Frequency  (rad/sec)’) 

axis([0,2,-2,3]) 

figure(gcf+l) 

plot(freq*2*pi,logl  0(zabcl  1)),  grid 

xlabel(‘Frequency(rad/sec)’) 

axis([0,2,-3,2]) 


Hmat.m 


%  This  program  is  used  for  the  experimental  data  obtained  using  the  I-DEAS 

%  software.  It  takes  the  data  stored  in  beamb  and  creates  an  FRF  matrix  as  well  as 

%  an  inverse  FRF  matrix  for  any  boundary  condition.  The  user  inputs  the 

%  frequency  range  and  the  nodes  to  be  used  in  the  a-set. 

% 

%  Written  by  C.P.  DeGregory 

clear 

load  beamb  %  Loads  data  file  containing  complex  amplitudes  for  each  frequency  data 
%  point  in  the  spectrum, 
hll  =  Hbase(:,l); 

%  Plot  the  FRF  for  the  free-free  beam  (There  are  1601  spectral  lines:  one  every  1.25  Hz) 

plot(0: 1 .25 :2000,20*loagl0(abs(hl  1  ))),grid; 
axis([0,500,-40,50]); 

first  =  input(‘Enter  first  frequency:  ‘); 
last  =  input(‘Enter  last  frequency:  ‘); 

start  =  (first/1 .25)  +  1 ; 
finish  =  (last/1 .25)  +1 ; 

%  Initialize 

icnt  =  0; 

H  =  zeros(5*((fmish-start)+l),5); 
ffeqp  =  [first:  1 .25  :last] ; 

Zap  =  zeros(length(freqp),l); 

aset  =  input(‘ Enter  nodes  to  be  used  in  a-set  (vector  form):  ‘); 

%  Build  H  matrix 
for  i  =  (start)  :  (finish) 

H((5*icnt-4):5*icnt,:)  =  [Hbase(i,l:5);  Hbase(i,6:10);  Hbase(i,lL15);... 

Hbase(i,  16:20);  Hbase(i,21:25)); 

%  Build  Haa  matrix 

Ha((length(aset)*icnt)-length(aset)+l:  (length(aset)*icnt), :  )  =  H  ((5*... 
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icnt-1)  +aset,  aset); 

Za((length(aset)*icnt)-length(aset)+l:  (length(aset)*icnt), : )  =  inv(Ha((length. . . 
(aset)*icnt)-length(aset)+l:  (length(aset)*icnt), : ); 

%  1,1  Element  for  plotting 

Zap(icnt)  =  Za(((length(aset)*icnt)  -  length(aset))  +1,1); 


end 

plot(freqp,  20*logl0(abs(Zap))),grid 
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B.  SUPPORT  PROGRAMS 


Assemble2Beams.m 


%  This  program  assembles  the  global  mass  and  stiffness  matrices  of  the  two  beams 

%  defined  in  “Build2Beams”  by  combining  element  k  and  m  matrices  solved  in  the 

%  function  “fbeamkm 
% 

%  Loop  over  the  two  beams: 

for  icnt_beams  =1:2; 

k=zeros(ndof,ndof); 

m=zeros(ndof,ndof); 

%  Loop  over  the  number  of  elements  in  the  structure: 


for  elnum=l:num_elements; 

dofl  =2*connect(elnum,  1 )- 1 ; 
dof2=dofl+l; 

dof3=2*connect(elnum,2)-l ; 
dof4=dof3+l; 

%  ...  set  up  beamel  function  call: 

l_temp=element_length;  %  Using  fixed  element  lengths 

ei_temp=element_EI(elnum,icnt_beams); 

m_temp=element_mass(elnum,icnt_beams); 

[kbeam,mbeam]=fbeamkm(l_temp,ei_temp,m_temp); 

k(dofl  :dof2,dofl  :dof2)=k(dofl  :dof2,dofl  :dof2)+kbeam(l  :2,1 :2); 
k(dofl  :dof2,dof3  :dof4)=k(dofl  :dof2,dof3  :dof4)+kbeam(  1:2,3 :4); 
k(dof3:dof4,dofl  :dof2)=k(dof3  :dof4,dofl  :dof2)+kbeam(3:4,l  :2); 
k(dof3  :dof4,dof3  :dof4)=k(dof3  :dof4,dof3  :dof4)+kbeam(3 :4,3 :4); 

m(dofl  :dof2,dofl  :dof2)=m(dofl  :dof2,dofl  :dof2)+mbeam(l  :2,1 :2); 
m(dofl  :dof2,dof3:dof4)=m(dofl  :dof2,dof3:dof4)+mbeam(l  :2,3 :4); 
m(dof3  :dof4,dofl  :dof2)=m(dof3  :dof4,dofl  :dof2)+mbeam(3 :4, 1:2); 
m(dof3  :dof4,dof3  :dof4)=m(dof3  :dof4,dof3  :dof4)+mbeam(3 :4,3 :4); 

%  end  loop  over  the  number  of  elements: 
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end 


%  Reassign  k  and  m  to  new  variables  and  add  lumped  masses 

if  icnt_beams  =  1 ; 

ka  =  k;ma  =  m; 
elseif  icnt_beams  ==  2; 
kx  =  k;mx  =  m; 
end 

%  End  icnt_beams  loop: 
end 

clear  dofl  do£2  dof3  dof4  ltemp  ei_temp  m_temp  icnt_beams 
clear  k  m  kbeam  mbeam  elnum 
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BeamA_Prompt.m  &  BeamProperties.m 


%  This  program  prompts  the  user  for  Beam  A  data.  The  first  inputs  required  are  the 
%  nominal  physical  properties.  The  short  program  file  “BeamProperties”  is  offered 

%  to  the  user.  The  data  listed  in  “BeamProperties”  is  the  data  that  was  used  for  all 

%  beams  throughout  this  thesis. 


disp('  ');disp(' '); 

disp('  Enter  nominal  physical  properties  for  first  beam') 
props_file  =  'n'; 

props_file  =  input('  Run  "BeamProperties.m"  script  to  define  nominal  properties  ?... 
(y/n) » \'s'); 

if  props_file  —  'y'; 

BeamProperties; 

else; 

total_length  =  input('  Enter  the  total  beam  length  in  inches:  '); 
num_elements  =input('  Enter  the  number  of  elements:  '); 
nominal_EI  =  input('  Enter  the  nominal  El  value  (lbf/inA2):  '); 
nominal_area  =  input('  Enter  the  nominal  cross  section  area  (inA2):  '); 
nominal  density  =  input(’  Enter  the  nominal  weight  density  (lbf/inA3):  '); 

end; 


for  icnt_elements  =  l:num_elements; 

connect(icnt_elements,l:2)  =  [icnt_elements,icnt_elements+l]; 
end 


ndofi=length(connect)*2+2;  %  this  is  unrestrained  beam  DOF 

elementjength  =  total_length/num_elements; 
element_EI(:,l)  =  nominal_EI  *  ones(num_elements,l); 
element_area(:,l)  =  nominal_area  *  ones(num_elements.  1 ); 
element_density(:,l)  =  nominal_density  *  ones(num_elements,l); 
element_mass(:,l)  =  (element_density  .*  element_area  *  element_length)/386.4; 


%  End  Data  input  for  1st  beam:  Copy  data  for  2nd  beam: 

element_EI(:,2)  =  element_EI(:,l); 
element_area(:,2)  =  element_area(:,l); 
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element_density(:,2)  =  element_density( : ,  1 ) ; 
element_mass(:,2)  =  element_mass(:,l); 

%  End  BeamA_Prompt.m 


% 

% 


%  BeamProperties.m 

%  This  is  the  file  to  load  nominal  beam  data. 

total_length  =  60; 
num_elements  =10; 
nominal_EI  =  5 . 5  e5 ; 
nominal_area  =  .75; 
nominal_density  =  .283; 
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BeamXPrompt.m 


%  This  program  prompts  the  user  for  Beam  X  data.  Specifically,  the  user  is  asked  to 
%  choose  whether  to  input  element  El  modifications  or  to  have  them  chosen 
%  randomly.  If  chosen  randomly,  the  data  is  saved  in  the  binary  (.mat)  file 
%  “changedata”  and  is  then  offered  as  a  choice  for  the  next  run  of  this  code.  By 

%  doing  this,  the  option  exists  to  create  2  free-free  beams  with  unknown  damage  and 

%  2  a-set  constrained  beams  with  the  same  (but  still  unknown)  damage. 


dispC  ');disp(' '); 

dispC  Modify  nominal  physical  properties  for  second  beam') 


%  Adjust  El  values  for  second  beam: 

option  =  'u'; 
optionl  =  'n'; 
changeJEI  =  'n'; 

change_EI  =  input('  Modify  single/range  element  El  values  (y/n)?  ','s'); 
option  =  input('  User  defined,  random,  or  neither  (u,r,n)?  ','s'); 
optionl  =  input('  Read  previous  modification  values  (y/n)?  ','s'); 
if  optionl  =  y 

load  changedata 
change_EI  =  'n'; 

end 

while  (change_EI  ~=  'n')  &  (option  ==  'u') 


disp('  Enter  element  label(s)  for  El  modification') 
disp('  Use  MATLAB  vector  format>  1  3  5:7  9  ') 
new_lbls  =  input('  »  ',’s'); 

new_lbls  =  eval(['[',new_lbls, ']’]);  %  Converts  string  to  vector  of  labels 


disp('  Enter  El  change  for  element  range') 

EI_change  =  input('  Enter  percentage  El  change  (+/-  %) '); 
element_EI(new_lbls,2)  =  element_EI(new_lbls,2)+... 
(Exchange/ 100)  *  element_EI(new_lbls,2); 

change_EI  -  input('  Modify  another  element  El  value  (y/n)?  ','s'); 

end;  %  end  while 
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if  option  =  V 

numels  =  fix(3*rand+l); 
for  i=l  rnumels 

new_lbls  =  fix(10*rand+l); 
format  short 
EI_change  =  -rand/2; 

element_EI(new_lbls,2)  =  element_EI(new_lbls,2)+... 
EI_change*  element_EI(new_lbls,2) ; 
end 


save  changedata  numels  new_lbls  element_EI 

end 


%  End  BeamX_Prompt.m 
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BeamSensPrompt.m 

%  This  program  prompts  the  user  for  design  variable  element  label  defintitions.  It  is 
%  called  by  the  program  “BeamSensitivity.” 

max_num_dv  =  50; 

o/o  ***************  ei  j)Y  ELEMENT  LABEL  PROMPTING  ***************** 
define_EI_dv  =  input('  Calculate  element  El  sensitivities  (y/n)? 
if  define JEI_dv  =  'y'; 

all_elements  =  input('  Define  El  DV"s  for  all  elements?  (y/n)  Vs'); 
if  all_elements  =  'y'; 

for  icnt_elements  =  1  :num_elements; 

EI_el_lbl(icnt_elements,l) =  icnt_elements; 

end; 

define_EI_dv  =  'n'; 

else;  %  all_elements  ==  no,  user  selects  elements 

for  icnt_dv  =  1  :  max_num_dv; 

sprintf('  Enter  element  label(s)  for  El  DV  #  %3i',icnt_dv) 
disp('  Use  MATLAB  vector  format>  1  3  5:7  9  —  "d"  to  finish') 
new_lbls  =  input('  » ','s'); 
if  new_lbls  ~=  'd'; 

new_lbls  =  eval(['[',new_lbls, ']']);  %  Converts  string  to 

%  vector  of  labels 

EI_el_lbl(icnt_dv,l  :length(new_lbls))  =  new_lbls; 

else 

break 

end; 

end;  %  for  icnt_dv 
end;  %  if  all_elements 

mass_el_lbl  =  [  ]; 

end; 

o/o  *************  end  ELEMENT  LABEL  PROMPTING  ********************* 
%  End  script. 
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Boundary_Conditions.ni 


%  This  program  prompts  the  user  for  boundary  condition  information,  creates  a 
%  vector  of  DOF  (with  respect  to  the  unrestrained  structure)  and  then  extracts  the 

%  rows  and  columns  of  the  complementary  DOF.  The  program  defines  vector 

%  "free_dof_set"  containing  the  list  of  unrestrained  DOF.  The  boundary  conditions 
%  for  the  given  system  are  applied  in  this  script. 


if  exist('free_dof_set')=0;  %  Build  ffee_dof_set  vector 

disp('  Select  a  boundary  condition  set:') 
disp('  ( 1 )  Clamped-free') 
disp('  (2)  Clamped-Clamped') 
dispC  (3)  Pinned-Pinned') 
disp('  (4)  User-Defined') 
disp('  (5)  Free-Free') 

BC  Choice  =  input(' »  Enter  choice: '); 

if  BC_Choice  ==  1 ;  %  Clamped-free 

free_dof_set  =  [3:ndof]; 
restraints  witch  =  'y'; 

elseif  BC_Choice  =  2;  %  Clamped-Clamped 

free_dof_set  =  [3:ndof-2]; 
restraint_switch  =  'y'; 

elseif  BC_Choice  =  3;  %  Pinned-Pinned 

ffee_dof_set  =  [2:ndof-2  ndof]; 
restraints  witch  =  'y'; 

elseif  BC_Choice  =  4;  %  User-Defined 

icnt_dof  =  0; 
adddof  =  'y'; 
while  add_dof  ==  'y'; 

bc_node  =  input('  Node  number  for  restraint  ?  "0"  to  end: '); 
if  bc_node  ==  0; 
break 

end; 

bc_coord  =  input('  Translation  or  Rotation  ?  (t/r)  ','s'); 
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icnt_dof  =  icnt_dof  +  1; 
if  bc_coord  ==  't'; 

bc_DOF(icnt_dof)  =  2*  bcjnode  - 1; 
elseif  mass_coord  ==  V; 

bc_DOF(icnt_dof)  =  2*  bc_node; 
end;  %  End  if-elseif  block 

end;  %  End  while  add_dof 

bc_boolean  =  ones(ndof,l);  %  [1  1  1  ...  icnt_dof] 

bc_boolean(bc_DOF)  =  zeros(length(bc_DOF),l);  %  Put  zeros  in  restrained  dof 
all_dofs  =  [  1  :ndof] ;  %  List  of  all  dof 

free_dof_set  =  all_dofs.*bc_boolean';  %  Extract  free  dof 
restraint_switch  =  'y'; 

elseif  BC_Choice  =  5;  %  Free-free  beam 

free_dof_set  =  [1  :ndof]; 
restraint_switch  =  'n'; 

end;  %  End  if-elseif  choice  block _ 

end;  %  End  exist  block 

free_dof_set(free_dof_set==0)=[  ] ; 

ka  =  ka(free_dof_set,free_dof_set); 
ma  =  ma(free_dof_set,free__dof_set); 
kx  =  kx(free_dof_set,free_dof_set); 
mx  =  mx(free_dof_set,free_dof_set); 

%  ***************  END  BOUND ARY_CONDITION S . M  ************ 
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c. 


FUNCTIONS  CALLED  BY  MAIN  &  SUPPORT  PROGRAMS 


function  [kbeam,mbeam]  =  fbeamkm(l_temp,ei_temp,m_temp) 

%  This  function  computes  the  elemental  stiffness  and  mass  matrices  defined  using 

%  Galerkin’s  Method  and  Hermitian  shape  functions.  [Ref  7] 

krowl  =  [12  6*l_temp  -12  6*l_temp]; 

krow2  =  [6*l_temp  4*l_tempA2  -6*l_temp  2*l_tempA2]; 

krow3  =  [-12  -6*l_temp  12  -6*l_temp]; 

krow4  =  [6*l_temp  2*l_tempA2  -6*l_temp  4*l_tempA2]; 

kbeam  =  (ei_temp/l_tempA3)*[krowl;krow2;krow3;krow4]; 

mrowl  =  [156  22*l_temp  54  -13*l_temp]; 

mrow2  =  [22*l_temp  4*l_tempA2  13*l_temp  -3*l_tempA2]; 

mrow3  =  [54  13*l_temp  156  -22*l_temp]; 

mrow4  =  [-13*l_temp  -3*l_tempA2  -22*l_temp  4*l_tempA2]; 

mbeam  =  (m_temp/420)*  [mrowl  ;mrow2;mrow3;mrow4]; 


function  [phi_normal,orth]=massnormaI(phi,mass) 

%  This  function  mass  normalizes  a  modal  matrix. 

a=size(phi); 
nummodes=a(l  ,2); 
phi_normal=zeros(size(phi)); 

% 

for  ii=l  :nummodes; 

modalmass(ii)=phi( :  ,ii)'*mass*  phi(:  ,ii) ; 
if  modalmass(ii)~=0, 

phi_normal(:,ii)=(l/sqri(modalmass(ii)))*phi(:,ii); 

else 

phi_normal(:,ii)=phi(:,ii); 

end 

end 

% 

%  do  the  orthogonality  calculation: 

% 

orth=phi_normar*mass*phi_normal*  1 00; 
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function  [omega, phi]=fmodes(k,m) 

%  This  function  returns  a  vector  containing  eigenvalues  (rad/sec)2  and  a  matrix 
%  containing  the  mass  normalized  mode  shapes.  The  mode  information  is  sorted  by 
%  frequency  in  ascending  order. 

[v,d]=eig(m\k); 

[omega, index]=sort(diag(d)); 

[phi,orth]=fmassnormal(v(:,index),m); 

format  long 
disp('Freqs  in  Hz.:') 
disp(' - ~~~ - ') 

disp(sqrt(omega(  1  :min(length(omega),  1 2)))/2/pi) 
disp(' - ') 


function  [oset]  =  fOset_from_Aset(ndof,aset); 

%  This  function  determines  the  complementary  o-set  from  a  set  [1 : 1  :ndof]  and  the 
%  subset  a-set  =  [x  x  x  ...]. 

% 

%  ndof:  Total  number  of  DOF.  Set  is  labeled  "nset". 

%  aset:  Retained  DOF  (proper  subset  of  [  1 : 1  :ndof]) 

nset  =  [l:ndof]; 

for  icnt  =  1  :  length(aset); 

indices(icnt)  =  find(nset  ==  aset(icnt)); 

end 

bool  =  ones(size(nset)); 
bool(indices)  =  zeros(size(indices)); 
oset  =  nset(fmd(bool>0)); 
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function[k,m]=fspringmass(conn, mspec, fix) 

%  This  function  assembles  the  stiffness  [k]  and  mass  [m]  matrices  for  an  assemblage 
%  of  springs.  The  connectivity  is  established  in  the  array  'conn',  the  boundary 

%  conditions  (fixity)  are  specified  in  the  array  'fix'.  The  mass  is  specified  in  the 

%  vector  mspec. 

% 

%  For  example,  to  build  the  following  system: 

% 

%  m  2m  m 

%  |-.////--*~////~*==////==* 

%  k  k  2k 

% 

%  the  connectivity  array  is  as  follows: 

%  [12  ' 

%  2  3 

%  3  4 

%  3  4] 

%  The  connectivity  array  is  therefore  of  dimension: 

%  (number  of  springs)*(2). 

% 

%  The  fixity  array  for  the  above  system  is 
% 

%  fix  =  [0111] 

% 

%  The  fixity  array  is  therefore  of  dimension: 

%  (l)*(number  of  nodes). 

% 

%  The  mass  distribution  is  given  by  mspec: 

% 

%  mspec  =  [12  1] 

% 

%  The  mspec  array  is  of  dimension: 

%  (l)*number  of  dynamic  (free)  DOF 


kspring=[l  -1;-1  1]; 
ksize=length(fix); 
k=zeros(ksize,ksize); 
m=zeros(length(mspec)); 

%  assemble  stiffness  matrix: 

for  spring  =  1  :size(conn,l); 

index  1  =conn(spring,  1 ); 
index2=conn(spring,2); 
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k(index  1 ,  index  1  )=k(index  1 , index  1  )+kspring(  1,1); 
k(index  1  ,index2)=k(index  1  ,index2)+kspring(  1 ,2); 
k(index2, index  1  )=k(index2,indexl)+kspring(2, 1 ); 

k(index2,index2)=k(index2,index2)+kspring(2,2); 

end; 

%  zero  rows  and  columns  corresponding  to  boundary  conditions 

indices_to_keep  =  fmd(fix); 

if  length(indices_to_keep)  ~=  0; 

k  =  k(indices_to_keep,indices_to_keep); 
end 

%assemble  mass  matrix: 

m=zeros(size(k)); 
for  i=l:length(m); 

m(i,i)=mspec(i); 

end 
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